AQ-A065  155 


UNCLASSIFIED 


ENVIRONMENTAL  RESEARCH  AND  TECHNOLOGY  INC  CONCORD  MASS  F/G  4/1 
MICROMAVE-INFRARED  RETRIEVALS. (U) 

APR  78  H K BURKE*  R K CRANE*  M 6 FOULER  F19628-T6-C-0135 


■ 


AFGL-TR-79-0019 


Unclassified 


ssified 


SECURITY 


ATION  OF  THIS  "AGE  flWi«n  note  Entered) 


'PORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


W.  OOVT  ACCESSION  NO. 


•jyr  >p_erio 
depart* 


RIOO  COVERED 


— I E ^ 

7c  _ on  a 


ICRO WAVE -INFRARED  RETRIEVALS  # 


IT  NUMBER 


S.  CONI  R ACT  OR  GRANT  NUMBERfa) 


ADDRESS 

Environmental  Research  & Technology,  Inct 

696  Virginia  Road 

Concord,  Massachusetts  01742 


i-  MONITORING  AGENCY  NAME  & ADDRESS/"  dill  a rant  from  Controlling  Ottl ca) 


So.  OECL .AMI  FI  C ATI  ON/  DOWN  OR  AOIN  G 


16.  DISTRIBUTION  STATEMENT  (ol  thia  Raport) 

Approved  for  public  release;  distribution  unlimited. 


It.  SURRl  EMENTARY  NOTES 


It.  KEY  WORDS  ( Conthnta  on  rarataa  alda  II  nacaaaary  and  Idantlly  by  block  numbar) 


Radiative  transfer 
Multiple  scattering 
Microwave  radiometry 


Infrared  radiometry 
Simulation  of  cloud  effects 


).  Amir  ACT  (Continue  on  nwh  elde  II  notoooorr  and  Identity  by  Slack  number) 

^ A fast  efficient  solution  to  the  plane -parallel  scattering  atmosphere 
radiative  transfer  problem  was  developed  for  use  in  simulating  the  effects 
of  clouds  and  rain  on  infrared  and  microwave  passive  remote  sensing  sys- 
tems. The  solution  was  based  upon  the  use  of  a variational -iterative  pror- 
cedure.  Application  of  the  variational -iterative  procedure  to  simulating 
the  effects  of  low  stratus  and  thin  cirrus  clouds  revealed  that  for  the  window 

channel  (channel  7)  of  the  DMSP  infrared  sensor  system,  multiple  scattering' 


EDITION  OF  I NOV  ••  It  OMOLtTE 

T|U> 


1.  REClR'FNT**  CATALOG  NUMBER 


17.  DISTRIBUTION  STATEMENT  (ol  the  obotroet  entered  In  Block  20.  II  dltterenl  from  Report) 


IS.  SECURI 

Unclassified 


— —i 


» 

< 


I 


'V 


- .•r’-r 

T * 


TABLE  OF  CONTENTS 


f 

r 

( 


Page 


ABSTRACT  1 

1 . INTRODUCTION  6 

1.1  Objectives  6 

1.2  Summary  of  Results  7 

2.  BACKGROUND  11 

2.1  Required  Retrieval  Accuracy  11 

2.2  Clear  Sky  Retrievals  13 

2.3  Cloud  Effects  15 

3.  A SOLUTION  TO  THE  RADIATIVE  TRANSFER  EQUATION  WITH 

SCATTERING  19 

3.1  Variational-Iterative  Method  19 

3.2  Extension  to  Anisotropic  Scattering  25 

3.3  Results  of  Model  Calculations  29 

4.  INFRARED  RADIANCE  SIMULATIONS  WITH  CLOUDS  57 

4.1  Software  Modification  57 

4.2  Cloud  Simulations  57 

5.  PARTIAL  CLOUDINESS  IN  THE  FIELD-OF-VIEW  67 

6.  CONCLUSIONS  69 

REFERENCES  70 


3 


' 


l 

1 

< 


« 


LIST  OF  ILLUSTRATIONS 


Figures  Page 

1A  Source  Function  Profile  37 

IB  Mean  Intensity  (I)  and  Diffused  Flux  (F)  Derived 

from  Source  Function  in  Figure  1A  38 

2A  Source  Function  Profile  39 

2B  Mean  Intensity  (T)  and  Diffused  Flux  (F)  Derived 

from  Source  Function  in  Figure  2A  40 

3 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

A at  All  Eight  IR  Channels  60 

4 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

B at  All  Eight  IR  Channels  61 

5 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

C at  All  Eight  IR  Channels  62 

6 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

D at  All  Eight  IR  Channels  63 

7 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

E at  All  Eight  IR  Channels  64 

8 Brightness  Temperatures  and  Radiances  for  Cloud  Type 

F at  All  Eight  IR  Channels  65 


4 


LIST  OF  TABLES 


Tables  Page 

1 RMS  Retrieval  Errors  - Midlatitude  Spring  Profiles  9 

2 A Comparison  Between  Exact  and  Calculated  Values  of 

the  Source  Functions  at  the  Top  and  Bottom  of  the 
Atmosphere  31 

3 Comparison  Between  Exact  and  Calculated  Values  of 
Intensities  at  the  Top  and  Bottom  of  the  Atmosphere 

at  Different  Angles  p 32 

4 Cases  of  Anisotropic  Scattering  with  Phase  Function  34 

5 Cases  of  Anisotropic  Scattering  with  Phase  Function  35 

6 Cases  of  Anisotropic  Scattering  with  Phase  Function  36 

7 Accuracy  of  the  VI  Technique  43 

8 Accuracy  of  the  VI  Technique  44 

9 One  Layer  Tropical  Cloud  46 

10  Five  Layer  Tropical  Cloud  I 47 

11  Five  Layer  Tropical  Cloud  II  48 

12  Summary  of  Atmospheric  Conditions  50 

13  Optical  Properties  for  Cirrus  Cloud  Layers  at  Different 

Thicknesses  at  Different  Wavelengths  52 

14  Reflection  (r)  and  Transmission  (t)  through  Different 

Cirrus  Cloud  Layers  at  X = 0.7  p 53 

15  Reflection  (r)  and  Transmission  (t)  through  Different 

Cirrus  Cloud  Layers  at  X = 3 p 54 

16  Reflection  (r)  and  Transmission  (t)  through  Different 

Cirrus  Cloud  Layers  at  X = 10  p 55 

17  Summary  of  Cloud  Models  and  Their  Optical  Properties 

at  12  pm  58 


S 


1 . INTRODUCTION 


1.1  Objectives 

The  objective  of  this  contract  is  to  provide  an  optimum  procedure 

% 

for  the  retrieval  of  atmospheric  temperature  profiles  from  satellite 
infrared  and  microwave  radiometer  data  in  the  presence  of  clouds.  Toward 
this  end,  a fast  and  efficient  solution  to  the  radiative  transfer  equa- 
tion was  developed  for  the  rapid  simulation  of  cloud  effects. 

The  work  under  this  contract  was  organized  into  a series  of  tasks 
some  of  which  were  specified  in  the  original  contract  (1  November  1975 
to  30  June  1976)  and  the  remainder  in  the  modified  contract  (1  July  1976 
to  30  April  1977).  Specifically,  the  tasks  were: 

\ 

Initial  Contract  - 

• Determine  an  optimum  combination  of  microwave  and  infrared  , 

data  to  infer  atmospheric  temperature  profiles  using  the 

Environmental  Research  § Technology  (ERT)  originated  micro- 
wave-infrared sounder  procedure, 

• Evaluate  the  utility  of  the  variational -iterative  (VI) 
algorithm  for  the  operational  determination  of  the  environ- 
mental parameters  of  clouds  and  precipitation. 

• Process  and  analyze  Skylab  digital  data  and  photography  to  1 

determine  techniques  for  using  near-infrared  data  to  dis- 
criminate between  clouds  and  snow  cover. 

Modified  Contract  - 

• Perform  simulation  analysis  and  use  the  ERT  statistical 
inversion  algorithm  to  evaluate  errors  in  inversion 
results. 

• Investigate  methods  to  optimize  inversion  results  for  a 
cloudy  atmosphere. 

• Modify  and  extend  the  V-I  technique  to  simulate  radiative 
transfer  in  a cloudy  atmosphere. 
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• Develop  a geometrical  model  for  partial  beam  filling. 

• Develop  cloud  distributions  to  be  used  with  the  geometri- 
cal model. 

These  tasks  fall  into  two  categories:  the  search  for  an  optimum 
method  for  the  inversion  of  infrared  and  microwave  data  in  the  presence 
of  clouds,  and  the  processing  of  Skylab  data.  Work  on  the  latter  task 
category  consisted  of  computer  processing  of  digital  data  from  three 
Skylab-4  passes  and  the  preparation  of  supporting  photographic  data. 

The  required  data  were  delivered  in  March  1976  and  used  by  AFGL  for 

further  analysis  (Bunting,  et  al,  1977).  Work  on  the  former  task  i 

category  will  be  considered  in  this  report;  work  on  the  latter  will  not 
be  considered  further. 

1.2  Summary  of  Results 

Under  a previous  contract.  Environmental  Research  5 Technology  (ERT) 
prepared  a set  of  computer  programs  for  the  Air  Force  Geophysics  Labora- 
tory (AFGL)  to  simulate  infrared  radiances  and  microwave  brightness 
temperatures  and  to  invert  the  simulated  radiometer  measurements  to 
estimate  atmospheric  temperature  profiles  (Weichel,  1976).  Two  inver- 
sion procedures  were  provided,  the  Statistical  Method  of  Parameter 
Estimation  (as  developed  by  ERT)  and  the  Minimum  Information  Method 
(provided  by  the  Air  Force  Global  Weather  Central).  The  microwave 
brightness  temperature  simulation  programs  were  developed  by  ERT  and  con- 
tained facilities  for  including  the  effects  of  precipitation  and  clouds. 

The  infrared  radiance  simulation  programs  as  provided  by  the  Air  Force 
Global  Weather  Central  (AFGWC)  did  not  have  facilities  for  the  inclusion 
of  cloud  effects.  Under  this  contract,  the  infrared  simulation  package 
was  modified  to  simulate  cloud  effects  using  the  variational -iterative 
technique  for  the  solution  of  the  radiative  transfer  equation  for  a 
scattering  medium. 

The  variational-iterative  (VI)  method  for  the  solution  of  the  radi- 
ative transfer  equation  for  a plane-parallel  scattering  atmosphere  with 
azimuthal  symmetry  was  developed  by  Sze  (1976)  for  the  special  case  of 
isotropic  scattering.  This  program  was  modified  to  include  a higher 
order  approximation  to  the  phase  function  (anisotropic  scattering)  and 
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adapted  for  use  in  the  infrared  radiance  simulation  package  operating  on 
the  CDC-6600  computer  at  AFGL. 

The  infrared  simulation  package  was  also  updated  to  calculate  radi- 
ances for  all  the  DMSP  infrared  sounder  channels.  At  the  same  time,  the 
program  was  updated  to  conform  to  the  latest  (early  1976)  operational 
version  of  the  AFGWC  retrieval  software.  In  the  process  of  working  with 
retrievals  using  the  updated  Minimum  Information  Method  as  used  by  AFGWC, 
rather  large  retrieval  errors  were  discovered.  The  sources  of  the  dis- 
crepancies were  investigated.  In  this  process,  a number  of  errors  and 
shortcomings  in  the  numerical  solutions  and  programming  logic  of  the 
AFGWC  programs  were  discovered.  AFGWC  was  advised  of  these  problems  by 
the  contract  monitor. 

The  corrected  infrared  simulation  package  (not  including  the  modi- 
fications to  include  cloud  effects)  and  the  microwave  simulation  package 
were  used  to  statistically  investigate  infrared,  microwave,  and  infrared 
plus  microwave  retrievals  of  temperature  profiles  for  clear  sky  condi- 
tions and  microwave  retrievals  for  cloudy  conditions.  Clouds  generally 
have  a catastrophic  effect  on  temperature  retrievals  obtained  using 
infrared  radiances.  In  most  retreival  methods,  cloud  contamination  is 
detected  early  in  the  process  and  contaminated  data  are  not  used  for 
retrieval.  Infrared  simulations  including  cloud  effects  were  made  with 
the  modified  program  only  for  conditions  that  might  not  be  detected  by 
the  cloud  screening  algorithms,  cases  involving  thin  cirrus  clouds  and 
low  level  stratus  clouds. 

The  retrieval  methods  were  compared  using  simulated  radiance  and 
brightness  temperature  values  calculated  using  an  ensemble  of  midlatitude 
(Washington,  D.C.)  spring  radiosonde  profiles.  Inversions  were  obtained 
for  clear  and  cloudy  (microwave  only)  conditions  with  and  without  addi- 
tive instrument  (measurement)  noise  using  the  Statistical  Method  for 
Parameter  Estimation  (SMPE)  and,  for  infrared  only,  the  Minimum  Informa- 
tion Method  (MIM) . The  results  are  shown  in  Table  1.  The  root  mean 
square  (rms)  differences  (errors)  give  a measure  of  the  difference 
between  the  actual  and  estimated  temperature  at  5 pressure  levels.  The 
results  show  that  for  clear  sky  conditions,  the  microwave  only  and  the 
infrared  only  SMPE  retrievals  yield  identical  results.  The  Infrared  plus 
Microwave  SMPE  retrievals  are  not  significantly  better.  Under  cloudy 
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conditions,  the  microwave  retrievals  were  not  significantly  different. 
These  results  show  that  microwave  measurements  are  superior  to  the 
infrared  measurements  for  situations  where  the  inherently  larger  field 
l.  of  view  of  the  microwave  antenna  does  not  Significantly  effect  the 

measurements.  Finally,  no  significant  difference  was  obtained  between 
the  MIM  and  SMPE  retrievals  although  MIM  appeared  to  do  slightly  better 
at  the  surface  and  poorer  at  400  mb. 
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2 . BACKGROUND 

2.1  Required  Retrieval  Accuracy 

Infrared  radiance  measurements  have  been  available  from  meteorolo- 
gical satellites  for  a number  of  years.  DMSP  satellites  have  been  rou- 
tinely providing  AFGWC  with  data  for  operational  estimation  of  tempera- 
ture profiles.  The  National  Weather  Service  (NWS)  also  routinely  obtains 
satellite  soundings  from  the  National  Oceanic  and  Atmospheric  Administra- 
tion (NOAA)  satellites.  In  the  case  of  the  National  Weather  Service, 
the  satellite  derived  temperature  soundings  obtained  near  synoptic  update 
time  (within  several  hours  of  0000  and  1200  UT)  are  used  to  provide  input 
to  their  numerical  weather  forecasting  model  for  data  sparse  regions 
over  the  oceans  (regions  without  adequate  radiosonde  data) . Recent 
analysis  of  the  utility  of  satellite  data  by  the  NWS  suggests  that  their 
forecast  model  works  as  well  without  satellite  data  as  with  satellite 
data  (Tracton  and  McPherson,  1977).  They  found  that  they  had  sufficient 
data  for  ocean  areas  to  adequately  update  the  forecast  model  and  satel- 
lite data  are  not  required. 

At  the  present  time,  most  computer  forecast  models,  the  users  of 
satellite  data, expect  the  data  to  be  similar  in  information  content  to 
radiosonde  data.  The  satellite  is  viewed  as  replacing  the  radiosonde  in 
regions  where  radiosonde  data  are  not  acquired.  The  wind  information 
provided  by  the  radiosonde  is  sought  from  observations  of  cloud  motion; 
the  temperature  and  humidity  data  from  infrared  radiance  or  microwave 
brightness  temperature  measurements.  With  the  expectation  that  satellite 
data  are  to  replace  radiosonde  data,  the  yardstick  for  comparison  and 
error  analysis  is  the  radiosonde.  Unfortunately,  the  satellite  observa- 
tions do  not  have  the  information  content  of  radiosonde  measurements. 

Data  are  obtained  from  a limited  number  of  infrared  and/or  microwave 
radiometer  channels  and  the  radiometer  data  are  often  redundant.  The 
radiosonde  measurements  on  the  other  hand,  are  essentially  independent 
at  different  heights;  data  at  each  reported  pressure  level  provides  new 
information.  The  satellite  data  cannot  compete  in  terms  of  information 
content. 

Temperature  profiles  derived  from  satellite  data  retreivals  appear 
to  be  quite  similar  to  the  radiosonde  observations.  The  rms  differences 
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reported  in  Table  1 are  between  2 and  5°K  (less  than  2 percent) . The 
success  of  the  temperature  retrievals  is  derived  not  from  the  adequacy 
of  the  radiometers  but  from  the  regularity  of  the  atmosphere.  The  sta- 
tistical method  for  parameter  inversion  utilizes  a large  number  of  atmo- 
spheric profile  measurements  and  effectively  selects  (calculates  in  a 
least  squares  sense)  the  profile  that  most  likely  would  produce  the 
measured  radiances  (or  brightness  temperatures) . The  Minimum  Information 
Method  constructs  the  profile  that  produces  the  smallest  perturbation 
to  the  firsc  guess  profile  required  to  match  simulated  radiances  calcu- 
lated using  the  perturbed  temperature  profile  to  the  measured  radiances. 

In  either  case,  additional  information  about  the  state  of  the  atmosphere 
- either  the  forecast  or  the  most  likely  state  - is  used  to  obtain  a 
retrieval . 

The  numerical  forecast  models  also  produce  temperature  profile  esti- 
mates based  upon  prior  measurements  and  simplified  equations  of  f*luid 
dynamics.  The  current  state  of  the  art  in  numerical  forecasting  is  a 
temperature  profile  with  smaller  errors  (differences)  than  those  reported 
in  Table  1 (Smagorinsky,  1969).  At  any  synoptic  time,  the  forecast  temp- 
erature profile  has  smaller  rms  errors  than  the  satellite  retrievals.  The 
forecast  errors  are,  however,  correlated  and  a forecast  model  with  no 
temperature  profile  measurement  input  will  slowly  depart  from  the  actual 
state  of  the  atmosphere,  the  rate  of  departure  depending  upon  the  parti- 
cular model. 

At  the  present  time,  with  no  massive  restructuring  of  the  numerical 
forecast  models  and  their  objective  analysis  schemes,  satellite  data  are 
not  competitive  with  radiosonde  data  or  with  forecast  data  in  limited 
data  sparse  regions  over  limited  time  intervals.  As  a part  of  the  First 
Global  GARP  Experiment  (FGGE) , estimates  have  been  prepared  for  the  maxi- 
mum temperature  errors  allowable  for  use  in  numerical  forecasting.  The 
estimated  maximum  error  is  1°  rms  (Bengtsson,  1975),  a value  not  currently 
achievable  using  satellite  data  even  in  the  absence  of  cloud  contamina- 
tion. Satellite  data  may  still  be  of  importance  for  military  application 
for  instances  in  which  no  radiosonde  data  are  available.  In  this  case, 
the  satellite  data  may  provide  the  only  connection  between  the  forecast 
model  and  reality.  Unfortunately,  through  the  strong  correlation  between 
the  first  guess  profile  and  the  retrieved  values,  the  forecast  model  may 
still  rapidly  depart  from  reality  when  MIM  is  used  to  provide  the 
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retrieval.  High  spatial  resolution  satellite  data  may  also  be  useful  in 
the  study  of  the  spatial  variation  of  integrated  properties  of  the  atmo- 
sphere such  as  total  precipitable  water.  This  aspect  of  the  utility  of 
satellite  data  is  not  considered  in  this  report. 

Modified  objective  analysis  schemes  may,  in  the  end,  be  required 
before  satellite  data  are  useful  for  input  to  numerical  models.  Although 
many  possible  atmospheric  profiles  may  produce  the  observed  radiance 
values,  a much  smaller  number  of  profiles  are  consistant  both  with  the 
spatially  developing  meteorological  fields  and  with  the  observed  hori- 
zontal variation  in  radiance.  The  further  reduction  of  the  number  of 
possible  profiles  will  reduce  the  rms  differences  between  estimated  and 
actual  temperature  values.  With  a statistical  parameter  estimation  pro- 
cedure operating  on  entire  radiance  maps  and  forecast  fields,  the  resul- 
tant temperature  estimation  error  may  drop  below  the  value  (or  any  value) 
suggested  as  a prerequisite  for  the  use  of  satellite  data.  With  the 
continued  development  of  statistical  data  assimilation  techniques  for 
objective  analysis  (Bengtsson,  1975),  opportunities  may  arise  for  impro- 
ving the  performance  of  temperature  data  retrieval  schemes  utilizing 
satellite  measurements. 

2.2  Clear  Sky  Retrievals 

The  programs  required  to  provide  clear  sky  retrievals  using  simu- 
lated data  were  provided  to  AFGL  by  ERT  under  a previous  contract.  The 
simulation  routines  and  retrieval  methods  are  described  by  Weichel  (1976). 

The  infrared  simulation  routine  was  modified  and  updated  under  this  con- 
tract to  conform  with  the  latest  version  (early  1976)  of  the  sane  routine 
used  operationally  in  the  MIM  retrieval  package  at  AFGWC.  The  program 
was  also  used,  at  the  request  of  the  contract  monitor,  to  obtain  retrie- 
vals from  a set  of  DMSP  measurements  rather  than  from  simulated  radiance 
values.  In  the  course  of  trying  to  match  the  observed  and  retrieved 
temperature  values,  it  was  discovered  that  better  results  could  be 
obtained  by  slightly  modifying  the  transmittance  values.  These  results 
suggest  that  the  filter  characteristics  used  to  prepare  the  transmittance 
values  for  the  simulation  package  may  be  in  error. 

The  results  obtained  from  the  infrared  simulation  package  are  for  a 
hypothetical  infrared  radiometer  system  that  has  the  filter  character- 
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istics  used  to  generate  the  transmittances  for  simulation.  To  the  extent 
that  the  simulation  package  represents  an  approximation  to  the  complex 
infrared  absorption  and  scattering  processes  that  occur  in  the  atmosphere, 
the  simulation  results  are  for  the  hypothetical  system  and  atmosphere. 
Difficulties  with  the  simulation  of  radiances  obtained  using  NOAA  satel- 
lites has  caused  the  National  Environmental  Satellite  Service  (NESS)  to 
abandon  the  MIM  approach  to  retrieval,  which  requires  a simulation  of 
the  atmosphere  and  instrument  response,  in  favor  of  the  SMPE  approach 
using  measured,  not  simulated,  radiances  to  construct  the  D matrix  (see 
Weichel,  1976  for  a definition  of  terms;  see  Smith  and  Woolf,  1976  for  a 
discussion  of  the  NESS  retrieval  system  for  Nimbus-6) 

An  investigation  of  the  errors  that  result  from  using  the  SPME 
retrieval  technique  and  combinations  of  infrared  and  microwave  measure- 
ments was  made  to  determine  whether  the  microwave  plus  infrared  measure- 
ments provide  more  information  than  either  data  type  alone.  A midlati- 
tude spring  data  set  was  used  for  the  simulations  and  temperature  pro- 
file retrievals.  Radiosonde  data  for  March,  April,  and  May  for  the 
years  1958  through  1962  were  used  for  the  analysis.  A detailed  discus- 
sion of  the  data  set  and  the  method  used  to  extend  the  data  set  to  reach 
to  a 10  mb  pressure  level  is  given  in  Weichel  (1976) . 

The  simulations  were  made  for  an  over  water  or  ocean  situation.  The 
ocean  surface  was  assumed  to  be  a black  body  at  infrared  frequencies 
(emissivity  equal  to  1.).  At  microwave  frequencies,  a more  complex  ocean 
surface  model  was  used  that  included  both  salinity  and  sea  surface  rough- 
ness effects  (Gaut,  et  al,  1972).  For  both  data  sets,  the  sea  surface 
temperature  was  assumed  to  have  a random  variation  about  the  surface 
(1000  mb)  level  temperature  value.  The  sea  surface  temperature  deviation 
was  normally  distributed  with  a 1.5°K  standard  deviation.  For  the  micro- 
wave  simulations  the  wind  speed  was  assumed  to  have  a 3.8  m/s  mean  value 
and  a 2.1  m/s  standard  deviation.  The  salinity  was  assumed  fixed  at  0.66 
moles/liter. 

The  results  of  the  simulations  with  and  without  additive  radiometer 
measurement  noise  are  listed  in  Table  1.  The  infrared  radiometer  data 
included  observations  (simulated)  from  8 channels.  The  microwave 
sounder  included  observations  (simulated)  from  7 channels.  The  frequen- 
cies and  weighting  functions  for  each  channel  are  specified  in  Weichel 
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(1976).  The  channel  characteristics  are  for  the  DMSP  infrared  sounder 
SSE  sensor  package  plus  the  microwave  sounder  package  to  be  used  on 
future  DMSP  satellites.  The  results  are  nearly  identical  to  those 
expected  for  the  SSH  package  when  the  535.35  cm  * F channel  is  included 
with  the  E channels. 

For  comparison.  Smith  and  Woolf  (1976)  reported  a similar  study  for 
the  17  channel  High-resolution  Infra-Red  Sounder  (HIRS)  and  5 channel 
SCAnning  Microwave  Spectrometer  (SCAMS)  on  the  Nimbus-6  satellite.  For 
clear  sky  conditions,  they  obtained  results  similar  to  those  presented 
in  Table  1,  simulated  errors  between  1°  and  2°K  at  altitudes  above  800  mb 
for  7 infrared  channels  near  15  pm  and  for  the  5 microwave  channels  plus 
the  same  7 infrared  channels.  They  did  not  analyze  microwave  only  data. 
Their  results  however  differ  from  ours  near  the  surface  where  their  in- 
version errors  are  all  less  than  0.5°K.  They  were  not  explicit  about  the 
method  used  to  simulate  surface  temperature  or  on  the  role  played  by 
surface  temperature  in  the  inversion  process.  At  one  point  they  comment 
that  surface  temperature  may  be  treated  as  a known  providing  an  extra 
channel  of  information.  They  therefore  may  have  artificially  constrained 
the  surface  temperature  values. 

In  summary,  for  the  DMSP  satellites  and  for  Nimbus-6,  the  cloud-free 
retrieval  errors  are  of  the  order  of  2°K  except  at  the  surface.  Smith 
and  Woolf  found  little  difference  between  using  infrared  alone  and  using 
infrared  plus  microwave  data.  We  find  little  difference  between  these 
two  cases  and  between  them  and  using  microwave  data  alone  whether  the 

f 

data  are  noise-free  or  simulated  as  contaminated  by  a reasonable  level 
of  instrument  noise. 

2.3  Cloud  Effects 

Clouds  may  severely  affect  infrared  observations  while  causing 
little  effect  at  microwave  frequencies.  Cloud  contamination  is  a major 
problem  which  affects  all  infrared  only  inversion  schemes.  Cloud  con- 
tamination detection  algorithms  are  used  to  remove  the  effects  of  clouds 
either  by  rejecting  contaminated  data  or,  in  the  case  of  NESS  retrievals 
(Smith  and  Woolf,  1976)  by  estimating  equivalent  cloud-free  radiances. 

Cloud  contamination  models  have  generally  been  straight  forward. 

At  infrared  wavelengths,  the  clouds  are  assumed  to  have  emissivities  near 
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one  and  large  optical  depths.  For  completely  overcast  conditions,  the 
clouds  provide  a lower  "surface"  and  observations  may  be  made  down  to 
cloud  level.  For  partly  cloudy  conditions,  the  clouds  obscure  only  part 
of  the  surface.  The  radiances  from  the  clouds  combine  with  those  from 
the  surface  in  proportion  to  their  areas.  In  effect,  two  radiative 
transfer  problems  are  solved,  one  for  absorption  by  the  atmosphere  above 
the  surface  and  the  other  for  absorption  by  the  atmosphere  above  a cloud. 
Each  are  handled  as  plane-parallel  problems  with  a lower  surface  having 
unity  emissivity  (see  Crane,  1976  for  notation). 

ll  = Bv[T(Ps)1VPs)  ' /PSMT(p)l  ^-dP  (D 

o “ 

:v  = Bv[T(pc)l0v(Pc)  ' /PcBvtTCP)]  (2) 

o v 

s 

where  I is  the  radiance  at  the  satellite  for  the  lower  surface 

v 

corresponding  to  the  earth's  surface,  s refers  to  surface 

c 

I is  the  radiance  at  the  satellite  for  the  lower  surface 

v 

corresponding  to  cloud  top,  c refers  to  cloud  top 

BV[T(P)]  is  the  Planck  radiance  (see  eq.  63) 

T(P)  is  temperature  (kinetic) 

P is  pressure 

v is  frequency 

6(p  ) = e Tv  = transmittance 

s 

s 

t is  optical  depth  = / y (y)dy 

o 

y (y)  is  the  absorption  coefficient  at  y along  the  path  from 
o to  s,  from  the  height  of  Ps  to  the  satellite. 


If  N represents  the  fractional  area  with  clouds  within  the  field  of  view 
of  the  radiometer,  then 


From  these  equations,  it  is  evident  that  the  radiance  values  will 
be  lower  if  clouds  appear  within  the  field  of  view.  Cloud  contamination 
may  be  detected  by  assuming  that  adjacent  scan  positions  for  a scanning 
radiometer  observe  different  areas  having  different  amounts  of  cloudiness 
but  having  identical  I®  radiances.  If  one  scan  position  has  no  cloud 
contamination,  it  will  have  the  highest  radiance  value.  Using  the  window 
channel  which  primarily  measures  surface  temperature  (and  water  vapor) 
the  effect  of  clouds  will  be  the  largest  due  to  the  greatest  change  of 
radiance  values  between  viewing  the  surface  and  viewing  clouds.  The 
cloud  detection  algorithm  that  is  used  by  AFGWC  assumes  that  the  highest 
radiance  value  for  a scan  area  is  for  the  cloud-free  location  and  the 
cloud-free  radiances  for  the  scan  area  are  taken  from  the  assumed  cloud- 
free  location.  This  technique  will  not  detect  and  correct  for  cloud 
contamination  if  all  scan  positions  view  areas  with  some  clouds. 

NESS  uses  the  data  from  adjacent  scan  spots  to  compute  the  frac- 
tional area  with  clouds  and,  from  the  computed  N values,  calculate  the 
equivalent  cloud  free  radiances.  This  can  only  be  done  if  two  additional 
pieces  of  information  are  known,  the  radiances  in  one  of  the  channels 
for  cloud-free  conditions  and  the  relative  heights  of  clouds  in  adjacent 
scan  positions.  If  the  sea  surface  temperature  is  known  or  adequately 
modeled  for  the  time  of  the  observations,  then  the  radiance  value  for 
the  window  channel  may  be  estimated  (using  climatological  data  to  account 
for  atmospheric  effects  such  as  the  additional  absorption  and  emission  by 

water  vapor),.  If,  further,  the  heights  of  the  clouds  in  adjacent  scan 

j c s 

positions  are  assumed  to  be  identical,  then  the  I valves  and  I values 
r v v 

for  adjacent  scan  spots  are  known.  Using  equation  (3),  the  ratios  of 
cloud-free  areas  for  adjacent  scan  spots  can  be  calculated  and  from  that 
the  equivalent  clear  sky  radiance  values.  Letting  N*  = the  ratio 

of  cloud-free  areas  (Smith  and  Woolf,  1976)  then 


where  I , and  I „ are  the  observed  radiances  for  the  adjacent  scan 
v,l  v, 2 

spots.  Equation  (4)  is  used  to  solve  for  N*  for  the  one  channel  for 
which  I*  is  known  (calculated  using  the  assumed  surface  temperature  and 
radiative  transfer  model). 


- • - ■ - 


This  method  or  any  other  for  the  correction  of  infrared  radiance 
data  using  infrared  observations  alone  is  subject  to  a large  number  of 
uncertainties.  The  microwave  observations  provide  a means  for  removing 
the  uncertainty.  Smith  and  Woolf  use  the  microwave  data  in  their  ver- 
sion of  the  SMPE  to  estimate  N and,  in  turn,  to  estimate  cloud-free 
radiances.  From  the  simulation  results  reported  in  Table  1,  it  is  evi- 
dent that  microwave  data  alone  will  accomplish  the  same  results  with 
little  difference  in  the  rms  retrieval  error.  The  combined  data  may 
however  be  useful  for  estimating  cloud  parameters  such  as  cloud  height 
and  fractional  cloud  area  (assuming  equation  (2)).  Both  cloud  height 
and  fractional  area  are  additional  parameters  that  can  be  obtained  from 
SMPE  inversions  of  combined  infrared  and  microwave  data.  At  this  junc- 
ture, an  adequate  cloud  simulation  package  is  required  to  provide  radi- 
ance estimates  for  different  cloud  height  and  fractional  area  values. 

The  required  cloud  simulation  package  is  considered  in  the  next  three 
sections. 

The  use  of  combined  infrared  and  microwave  data  for  temperature  and 
cloud  parameter  retrievals  requires  matched  observing  areas  at  infrared 
and  microwave  frequencies.  The  microwave  systems  however  inherently 
observe  larger  areas  due  to  physical  limitations  on  antenna  beamwidths. 

If  the  observing  areas  at  infrared  and  microwave  frequencies  are  identi- 
cal, then  the  microwave  data  alone  are  adequate  for  temperature  retrieval. 
If  large  differences  exist  between  the  observing  areas  then  the  microwave 

data  are  not  useful  in  solving  the  cloud  contamination  problem  although 

s 

the  microwave  data  should  provide  better  estimates  of  I than  estimates 
from  surface  temperature  models. 


3.  A SOLUTION  TO  THE  RADIATIVE 
TRANSFER  EQUATION  WITH  SCATTERING 


The  inclusion  of  clouds  in  the  radiative  transfer  equation  compli- 
cates the  problem  due  to  the  addition  of  scattering  by  the  cloud  particles. 
Instead  of  the  nonlinear  quadrature  formula  represented  by  equations  (1) 
and  (2),  the  radiative  transfer  equation  becomes  an  integral  equation 
incorporating  the  effects  of  multiple  scattering.  Numerical  solutions 
to  the  multiple  scattering  equations  have  been  available  for  some  time 
but  the  procedures  have  been  time  consuming  and.  costly  to  operate.  For 
example,  Kattawar  and  Plass  (1969)  presented  the  results  of  a Monte  Carlo 
simulation  of  the  multiple  scattering  from  water  clouds  at  visible  and 
near  infrared  wavelengths.  Sze  (1976)  recently  developed  a rapid  method 
for  solving  the  multiple  scattering  equation  for  isotropic  scattering  in 
a plane-parallel  medium  (the  scattering  properties  vary  only  with  height 
above  a plane  earth).  With  the  variational -iterative  method  developed 
by  Sze,  computer  times  are  short  enough  and  costs  low  enough  to  permit 
the  inclusion  of  multiple  scattering  effects  in  the  simulation  of  infrared 
radiances.  The  details  of  the  variational -iterative  solution  as  developed 
by  Sze  and  extended  under  this  contract  to  include  anisotropic  scatter- 
ing are  discussed  in  this  section.  The  software  package  provided  to 
AFGL  and  the  results  of  simulations  will  be  discussed  in  the  next  section. 


3.1  Variational -Iterative  Method 

The  radiative  transfer  equation  for  a plane-parallel  scattering 
atmosphere  is  given  by  (Chandrasekhar,  1960) : 

v = I v (t , p , d>)  - Jv(t,w,$) 


. +1  2ir 

where  Jv(t,w,<J>)  = / / p (w , ; p 1 , «(> 1 ) I (t , u 1 , <|> 1 ) dw  1d<J> 1 

-1  o p 

♦ Jq(t,V,<P) 

u is  the  cosine  of  the  zenith  angle 
4 is  the  azimuth  angle 
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J is  the  source  function 


Jq  is  the  primary  exitation  source  function 
[ p is  the  phase  function 

and  the  remaining  symbols  were  defined  above.  For  a scattering  plane- 
parallel  atmosphere  with  the  primary  source  emission  from  atmospheric 
gasses  or  for  a plane-parallel  atmosphere  with  an  external  source  (the 
sun  say)  and  isotropic  scattering,  the  problem  has  axial  symmetry  about 
the  zenith  direction  and  the  dependence  on  azimuth  can  be  removed: 

= I(t,p)  - J(T,p)  (5) 

J(t,p)  = J0(x,u)  + j / pOi.P^Ux.p^dp1  (6) 

where  the  explicit  dependence  on  v has  been  suppressed. 

The  phase  function  is  normalized  to  represent  the  total  energy 
scattered  by  a single  scatter  relative  to  the  sum  of  the  energy  absorbed 
and  scattered  by  the  particle: 

jr  / / p(w,«fr;Pl»'l>1)dw1d«|>1  = o)  < 1 (7) 

"♦TT  -J  O ° 

• 

where  is  the  single  scattering  albedo.  For  an  isotropic  scatterer, 
p = constant  = The  primary  exitation  source  function  for  atmospheric 

thermal  emission  is 

Jq(t)  - (1  - wo(t))B[T(t)]  (8) 

where  B[T(t)]  is  By,  the  Planck  function.  In  this  representation  of  the 
radiative  transfer  equation,  the  distance  from  the  top  or  bottom  of  the 
atmosphere  is  measured  in  units  of  optical  depth,  t,  rather  than  distance 
or  pressure.  The  point  properties  of  the  medium  such  as  the  single 
scattering  albedo  are  then  functions  of  t. 
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tCO  = / Ye(C)d C 
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ID  (O  = ~ 

o'-  YeU) 


(10) 


where  Yg  is  the  volume  extinction  coefficient  (cross  section  per  unit 
volume) 

Y is  the  absorption  coefficient  as  defined  above  and  C = height. 

2L 

Equation  (5)  reduces  to  (1)  above  when  Yg  = Ya  or  u = 0. 


From  equations  (7),  (8)  and  (9),  the  source  function  can  be  written 


as 


f N T * 

J(T)  = 1 - <do(t))B(T(t))  + I J(T)E1(|r-T*|)dT 

♦ 22£! 


(ii) 


where  E^  is  the  exponential  integral  of  the  first  order  and 

t*  is  the  total  optical  depth  such  that  I(t*,u)  is  the  outgoing 
(upwelling)  intensity  at  the  lower  boundary  (surface). 

The  exponential  integral  En  of  the  nth  order  is  defined  as 

En(x)  = /1e'x/pun'2du;  n = 1,  2-.. 
o 

I (x* , u)  has  two  contributions:  (1)  surface  emission  and  (2)  surface 
reflection.  For  a surface  reflectivity,  R,  and  temperature,  T , 


IjCt-.p)  = (1  - R)B(Ts)  . 


(12) 


I2(T*,lO  = 2*R*  / I(T*,-U1)w1dp1  (13) 

o 

where  I(T*,-p!)  is  the  downward  intensity  at  the  surface  in  the  direc- 
tion corresponding  to  u1. 


i 
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Equation  (13)  is  based  on  the  assumption  of  a Lambertian  surface;  there 
will  be  an  "isotropic"  reflection  independent  of  the  incident  angle  on 
the  surface.  Equation  (13)  can  further  be  written  as 


= t)  MT*-  ^dt  0*) 

where  E 2 is  the  second  order  exponential  integral.  Combining  (11),  (12) 
and  (13)  we  obtain  the  expression  for  the  source  function: 

J(T)  = (1  - <oo(t))B(T(t))  ♦ I WEjdt-xlJdt 

° (15) 

♦ * E2(T*-r)[(l-R)B(Tgr)  + 2-R  • /TJ(t)E2 (r*-t)dt] 


The  outgoing  intensity  at  the  top  of  the  atmosphere  then  can  be  expressed 


T * 

1(0,  v)  = / J(t)e‘t/wdt/u  + (l-R)B(T  )e'T*/lJ 
0 T* 

+ 2-R-[/  J(t)E  (T*-t)dt]e'T*/p 


The  three  terms  on  the  right  hand  side  of  equation  (16)  are 

(1)  the  upward  emission  from  the  atmosphere  (including  clouds), 

(2)  the  emission  of  the  attenuated  background  surface,  and 

(3)  the  reflection  of  the  downward  atmospheric  emission  from 
the  surface. 

The  variational -iterative  (VI)  approach  (Sze,  1976)  is  used  to 
solve  this  system  of  equations.  The  variational  method  depends  on  find- 
ing the  "extremum"  of  a certain  functional;  an  a priori  form  is  used 
for  the  unknown  function  and  the  coefficients  are  found  from  a set  of 
minimizing  conditions.  In  essence,  this  method  provides  a direct  way 
for  constructing  an  approximate  solution  for  the  source  function.  The 
atmosphere  is  divided  into  subintervals  and  the  source  function  is  approx- 
imated as  a combination  of  step  functions  in  different  intervals.  The 
advantages  of  this  technique  are  that:  (1)  it  is  fast  and  requires 
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little  computational  time  to  achieve  satisfactory  accuracy,  and  (2)  it 
allows  vertical  inhomogeneity  and  the  inclusion  of  surface  reflection. 

The  VI  technique  provides  a direct  method  for  constructing  an 
approximate  solution  to  the  integral  equation  (15)  for  the  source  func- 
tion. An  approximate  source  function  can  be  expressed  as 


Ja(T)  = Ua(T)\|u.o(T) 


where 


UaW  * l CiVi(,) 
1=1 


and  the  V. (t)  are  known  trial  functions.  The  choice  of  trial  functions 
i 

plays  an  important  role  in  the  ultimate  success  of  the  variational 
method  (Kourganoff,  1963).  In  the  variational  solution  employed  by  Sze, 
simple  step  functions  were  chosen  as  the  trial  functions.  This  choice 
(1)  makes  it  simple  to  perform  the  integrals  required  in  (15)  and  (2) 
the  intervals  can  be  chosen  to  resemble  multiple  cloud  layers,  the 
weights  \fuT”  for  each  layer  represent  the  average  source  function  in 
that  layer  where  ok  is  the  single  scattering  albedo  for  the  layer. 

The  total  optical  depth  t*  is  divided  into  N-l  intervals  with  ok 
constant  over  each  interval.  The  trial  function  then  is 


'i«-C 


1 T . < T < T . . 
3 ~ J+l 

0 otherwise 


The  then  are  solutions  of  the  algebraic  equation 


where 


ll 

y m. .c.  = f. 

j=i  ^ •*  1 
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D.(t)  = / ^EpT-t^dt 

J T 


,1-wi 


fj  = (-^-)B(T(Ti))  + -y-  E2(T*-Ti)  (1-R) B(TS) 
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The  variational  solution  is  an  approximation  to  the  actual  solution 
which  is  correct  at  least  at  one  level  within  each  layer  (Sze,  1976).  A 
smoothed  approximation  for  the  source  function  then  can  be  constructed: 


o(t) 


Jj(t)  = (1-u>0(t))B(T(t))  + E2(t*-t)(1-R)B(Ts) 


u , -.N-l  

+ [ C.<Ju.  (t)  D.(t) 

i=l  J J J 


(22) 


N-l 


+ “>o(t)  E (t*-t)  R l C.Ju).  (t)  J ^ + 1E  (T*-t)dt 
j=l  J 3 t. 

Since  the  smoothed  approximation  is  a summation  over  layers  with  oscil- 
lating residual  errors,  it  provides  a reasonable  first  estimate  of  the 
true  source  function  (Burke  and  Sze,  1977) . Improved  accuracy  may  be 
obtained  by  further  iterations  of  the  integral  equation  for  the  source 
function  as: 

Jn+i  (T)  - (1-‘Mt))B(T(t))  + e2(t*-t)(1-R)B(Ts) 

♦ f E1(|x-t|)Jn(t)dt  (23) 

0 

+ wo(t)  E2(t*-t)  R JT*E2(T*-t)Jn(t)dt 


The  residue  of  the  nth  iteration  is  defined  as 

Pn  - ^n-1 1 


A = 
n 


(24) 


By  specifying  the  maximum  residue  allowed,  the  iteration  process  then 
brings  the  source  function  to  desired  accuracy.  The  number  of  iterations 
for  the  system  also  depends  on  the  choice  of  the  number  of  step  functions 
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in  the  matrix  formalism  (the  number  N,  see  eq.  17).  For  example,  we  con- 
sider a fully  scattering  (u)q  = 1)  medium  with  t*  = 3.5.  By  dividing  the 
medium  into  two  step  functions  (N=2)  it  takes  seven  iterations  to  bring 
the  maximum  residue  in  the  source  function  to  less  than  1%.  However,  by 
representing  the  medium  by  four  step  functions  it  takes  only  two  itera- 
tions to  obtain  the  same  accuracy. 

3.2  Extension  to  Anisotropic  Scattering 

Rain  and  cloud  particles  range  in  size  from  smaller  than  a wavelength 
to  many  wavelengths  at  microwave  and  infrared  frequencies.  Over  this 
size  range,  the  energy  radiated  by  scattering  is  not  scattered  isotropi- 
cally in  all  directions  but  is  scattered  selectively  into  forward  and 
backward  lobes.  The  relative  intensity  of  the  scattering  in  any  direc- 
tion is  described  by  the  phase  function  (equation  7).  In  general,  the 
phase  function  for  randomly  oriented  scatterers,  and  for  natural  or 
unpolarized  thermal  emission  can  be  expanded  in  a series  of  Legendre 
polynomials  as 

oo 

p(cos  0)  = £ a>£  Pf(cos  0)  (25) 

H=0 

where  P = Legendre  polynomial 
0 = scattering  angle. 

In  addition  to  the  simplest  case  of  isotropic  scattering,  two  phase  func- 
tions have  found  extensive  use  in  describing  particulate  scatter,  Rayleigh's 
phase  function 


p(cos  0)  = (1  + cos^0)  (2b) 

which  may  be  used  when  the  particles  are  small  compared  to  a wavelength 
and 

p(cos  0)  = 0)0(1  + x cos  0)  (-1  <_  x ± 1)  (27) 

which  is  often  used  as  a first,  2-term,  approximation  to  equation  (25) 
for  problems  with  particles  large  relative  to  a wavelength  having  signi- 
ficant differences  between  the  forward  and  backscattered  energies.  X is 
used  to  adjust  the  relative  magnitude  of  the  forward  and  backscattered 
energies.  For  applications  at  infrared  frequencies,  the  latter  approxi- 
mation was  used  to  characterize  anisotropic  scattering. 
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The  two  term  approximation  is  often  used  in  the  study  of  radiative 
transfer  in  planetary  atmospheres  because  it  iz  easy  to  handle  analyti- 
cally and  provides  a good  approximation  to  the  integrated  effects  caused 
by  predominantly  forward  and  backward  scattering  particles.  For  radia- 
tive transfer  problems  with  azimuthal  symmetiy  or  in  problems  employing 
azimuthally  averaged  radiances, 

,2  if 

J cos  0 d4> ' = p p- 
o 

hence  p(p,p')  = wo(l  + xfp') 


(28) 


Using  this  form  for  the  phase  function,  the  source  function  becomes 


J(T,p)  = [/  nT.^'idp*  + x(t)u/  I (t  ' ,p')p'dp'] 


-1 


-1 


(29) 


or 


where 


and 


+ J0(t,p) 


J(t,p)  = Jj  (t)  + p J2 (T) 


J1(T)  = /^(T.p'JdM*  ♦ Jq1(t) 


J2  (T  ) ' °9-  x(t)J  HT.p’Jp'dp'  + J02(t) 


-1 


(30) 


(31) 


(32) 


Employing  the  formal  solutions  of  (5)  for  I(t,p)  and  I(t,-p),  the 
outgoing  and  inward  intensities, 

I(t,v)  = JT  J(T,p)exp(-(t-T)/p)  dt/p 
T 

= JT  (Jj(t)  ♦ P J2(t))  exp  (- (t-x ) /p)  dt/p 

T 

and 

I(r,-p)  = / J(t,-p)  exp(-(t-t)/p)  dt/p 
o 

= j^fJ^t)  - P J2(t))  exp(-(T-t)/p)  dt/p 
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Preceding  in  a manner  similar  to  the  development  in  section  3.1,  the 
boundary  conditions  for  the  surface  are  incorporated  to  produce  a coupled 
set  of  integral  equations  for  J^(t)  and 

Jj(t)  = /T[E1  ( | T-t  | ) Jj(t)  + E2(  | T-t  | ) b (t-T)  J2(t)]dt 


♦ ( 1 -uo (t) ) B(T(t)  ) + (1-R)  B CTs)  E2  (t*-t) 


+ wo(t)  R E2(t*-t )/T  [E2 (r*-t) J i (t)  - E3(T*-t)  J2(t)]dt 


and  J2(t)  = x (r)/T* [E2 ( | T-t  | )b(t-T) Jj  ft)  + E3  ( | T-t  | ) J2 (t) ]dt 


0)  , A 

(l-H)B(Ts)  E3 (t*-t) 


+ wo(t)x(t)R  E3(T*-T)/T*[F.2(T*-t)J1(t)  - E3(T*-t)J2(t)]dt 


where  the  function  b(x)  is  defined  as 

,,  . ,- 1,  x<0 

b(x)  = { 1,  x>0 

In  matric  notation,  (33)  and  (34)  can  be  expressed  as 

/ji(x)\  mq(t)  jt*  Ei  ( | r-t  | ) E2(|T-t|)b(t-T)  /Ji(t)\ 

I.J2(t)J  2 o x(f)E2(|T-t|)b(t-T)  X (t)  E3  ( | T-t  | ) V J2(t)i  dt 


1-w0(T)  jB(T  ) + (1-R)  ( Ez(T*'T)  \ B(T  ) 

0S  2 Ix(t)E3(t‘-t) 


♦ u0(t).R  / E2 (t*_t)  ] /T*[E2(T*-t)Ji(t)  - E3(T*-t)J2(t)]dt 

\x(t)E3(t*-t)/  o 


The  variational-iterative  technique  is  used  to  obtain  a solution  to  the 
coupled  integral  equation  set  (35). 
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A similar  set  of  linear  algebraic  equations  can  be  written: 

2n 


y m..c.  = f. 

j=i  13  3 1 


(36) 


The  size  of  the  matrix  M is  2n  x 2n  where  n is  the  number  of  steps.  The 
matrix  can  be  divided  into  quadrants,  Ml,  M2,  M3  and  M4  such  that 


Ml.  . = M . 

1J  iJ 

M2. . = M. 

ij  i,n+j 


M3. . = M . . 
ij  n+1 , j 


M4 . . = M . 
i,j  n+i ,n+j 


where  i=l,  2,...n,  j=l , 2,....,n,  and 


Ml  . = 6.. (At  ) - + 1 /Tj  + 1 Ej  (|T-t |)dt  dx 

1 > J J L 


Ti  Ti 


M2.  . = 
i»J 


fi+1JT3+1  E2(|x-t|)b(t-x)dt  dx 


Ti  Ti 


. . = - x(Ti )/T1  + 1/T^  + 1 E2  ( | T-t  | )b(t-x)dt  dx 

i i = 6ii  CAT.)  - iilallii  x(^r^rj  + 1 E3(|x-t|)dt  dx 
1,3  13  3 z 


M4 


,ith  sij  ■ <o’,  i ; s 

The  forcing  terms  f^  are  defined  as 


Ax.  = x.  , - x. 
l l+l  l 


(37) 

(38) 

(39) 

(40) 


fj  = JTi+1[(l-a>0(T))B(Ts)  + (l-R)E2(x*-x)B(Ts)]dx 

for  i=l , 2, . . .n  , 


Ti 


(41) 


and  f.  = r 1 l£flIlIx(T)  (l-R)E3(x*-x)B(Ts)dx 


(42) 


for  i=n+l, . . . ,2n 
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3.3  Results  of  Model  Calculations 


The  VI  technique  was  used  to  prepare  results  that  could  be  compared 
with  previous  calculations.  This  was  done  (1)  to  test  the  VI  technique 
and  the  programs  developed  under  this  contract  and  (2)  to  provide  insight 
into  the  effects  of  multiple  scattering  on  atmospheric  emission  estimates. 

3.3.1  Solar  Diffusion  Results 

A considerable  literature  exists  on  the  solution  to  the  radiative 
transfer  equation  where  the  primary  exitation  term  describes  radiant 
energy  from  the  sun.  In  this  case  (see  equation  (8)  for  comparison) 

J0(t,p)  = p(p,4V1°')  exp(-T/u0) 


where  pq  is  the  cosine  of  the  incidence  angle  (zenith)  for  the  solar 
radiation.  For  anisotropic  scattering  an  azimuthally  averaged  intensity 
may  be  determined  using  the  axial  symmetric  formulation  (eq.  6)  and  the 
phase  function  approximation  introduced  in  equation  (28), 


J0(t,p)  = exp(-T/u0)  (I-x(t)uuo) 


(43) 


The  coupled  integral  equations  then  are 


Ei  ( | T-t | ) 
xE2(|x-t|)b(t-x) 


E2 ( I x-t | )b(t- 

XE3 ( | T-t | ) 


Wall) 

4 


exp(-x/p0) 


dt 

(44) 


The  coupled  equations  are  solved  using  the  VI  technique  as  developed 
above.  Once  J^fx)  and  ^(x)  are  determined,  the  azimuth  and  zenith  angle 
averaged  mean  intensity  (I(x))  and  vertical  flux  (irF(x))  throughout  the 
medium  are  immediately  obtained  as 

T(x)  = 1/2  J1  I(x,p)  dp  = -4-v  fJl  (t)  - exp(-x/y0)]  (45) 

-1  wo tT/  H 

,1 

and  F(x)  * 2 / I(x,u)  pd*' 

-1  (46) 

* [«^2 ("0  ♦ x(T)n0exp(-x/wo)] 
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X ■ 0 represents  the  special  case  of  isotropic  scattering.  In  this 
case  the  source  function  is  reduced  to  J(x,p)  = J^x),  independent  of  p. 
Therefore,  by  setting  X = 0 in  solving  equation  (44),  we  get  solutions 
for  the  isotropic  system.  The  results  are  compared  with  the  exact  solu- 
tions obtained  from  Chandrasekhar's.  X,  Y functions  (1960)  by  setting 
the  maximum  residue  (A  max  %)  in  the  source  function  to  less  than  0.02%. 
The  relations  of  X,  Y functions  to  intensity  and  source  function  are: 


I(o, p)  = 1/4  o)0  F — J2—  [X(p)  X (p0)  - Y (p)  Y(p0)]  (47) 

M Mo 

I(T*,-P)  = 1/4  0)0  F y [Y (p)  X(P0)  - X(p)  Y (p0) ] (48) 

J(o)  = 1/4  uo  F X(po)  (49) 

J(t*)  = 1/4  u)Q  F Y(p0)  (50) 


The  tabulated  values  of  X,  Y functions  for  different  u)qs  and  x*s  were 
taken  from  Carlstedt  and  Mullikin  (1966). 

The  cases  presented  are  x*  * 1,  p » 1 for  u)Q  = .5  and  1.  Table  2 
is  the  comparison  of  the  source  functions  at  the  top  and  bottom  of  the 
atmosphere.  It  shows  that  for  <*>0  = .5,  the  calculated  values  are  close 
to  error-free;  they  agree  with  the  exact  values  at  up  to  four  figures 
after  the  decimal  point.  For  the  conservative  scattering  case  (u>o  - 1), 
the  error  is  still  within  0.2% 

Table  3 shows,  under  the  same  atmospheric  conditions  as  in  Table  2, 
the  comparison  of  intensities  at  the  top  and  bottom  of  the  atmosphere  at 
different  angles  (p).  The  results  indicate  that  for  uq  * .5,  the  errors 
in  intensities,  both  outward  and  inward,  are  generally  less  than  or  in 
the  neighborhood  of  0.1%.  For  a>o  * 1,  the  error  still  lies  within  0.3%. 

From  these  results,  we  are  convinced  that  the  variational-iterative 
scheme  is  a very  accurate  method  in  solving  radiative  transfer  problems. 
This  lays  an  important  foundation  for  the  anisotropic  cases  because 
"exact"  solutions  are  not  available  for  these  cases.  We  now  can  control 
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TABLE  2 

A COMPARISON  BETWEEN  EXACT  AND  CALCULATED  VALUES  OF  THE 
SOURCE  FUNCTIONS  AT  THE  TOP  AND  BOTTOM  OF  THE  ATMOSPHERE 
CASES  CONSIDERED  ARE  t*  = 1,  p = 1 for  u = .5  and  1 


(a) 

0 

J 

Exact 

Calculated 

Error  (%) 

l 

J(o) 

J(t*) 

0.439345 

0.241511 

0.438763 

0.241123 

0.13% 

0.16% 

.5 

J(o) 

J(T*) 

0. 15327] 

0.065926 

0. 153242 

0.065924 

0.02% 

0.00% 

TABLE  3 


COMPARISON  BETWEEN  EXACT  AND  CALCULATED  VALUES  OF  INTENSITIES  AT  THE  TOP  AND 
BOTTOM  OF  THE  ATMOSPHERE  AT  DIFFERENT  ANGLES  p.  CASES  CONSIDERED  ARE 
t*  = 1,  p = 1 for  o)q  = .5  and  1 


fa) 

0 

I(t,  p) 

Exact 

Calculated 

Error 

1 

1(0,  1) 

0.269393 

0.268759 

0.23% 

1(0,  0.8) 

0.30680 

0.306064 

0.24% 

1(0,  0.6) 

0.353534 

0.352666 

0.24% 

1(0,  0.4) 

0.409034 

0.407976 

0.26% 

CM 

© 

* 

o 
> ' 

1 

0.458037 

0.456708 

0.29% 

I(t*,  -0.8) 

0.274833 

0.274301 

0.19%  * 

I(t*.  -0.6) 

0.306037 

0.305458 

0.19% 

I(t*,  -0.4) 

0.332541 

0.331954 

©\° 

oo 

rH 

d 

I(T*,  -0.2) 

0.325061 

0.324623 

0.13% 

o.s 

1(0,  I) 

0.076583 

0.076521 

0.08% 

1(0,  0.8) 

0.087701 

0.087625 

0.09% 

1(0,  0.6) 

0.101983 

0.101888 

0.09% 

1(0,  0.4) 

0-120045 

0.119913 

0.11% 

1(0,  0.2) 

0.140087 

0.139865 

0.16% 

1 (t * , -0.8) 

0.074490 

0.074478 

0.02% 

I (t*,  -0.6) 

0.082344 

0.082339 

0.01% 

/■—v 

H 

* 

1 

O 

0.088386 

0.088396 

0.01% 

I(J*,  -0.2) 

0.084772 

0.084824 

0.01% 

the  desired  accuracy  by  iterating  the  solution  until  the  residue  in  the 
source  function  is  small  enough  to  be  acceptable. 

A wide  range  of  different  atmospheric  conditions  (both  on  x*  and 
and  on  phase  function  p(y,iO  = u)Q(l  + XuvO  were  tested.  Shown  in 
Tables  4,  5 and  6 is  a summary  of  the  representative  cases.  The  optical 
depths  (t*)  chosen  are  0.25,  1 and  4,  with  single  scattering  albedos 
a>o  = 1 and  0.5.  Three  x's,  -1,  0.5  and  1,  were  used.  The  exact  solutions 
were  obtained  by  the  variational  method  with  successive  iterations  such 
that  the  maximum  residue  in  the  source  function  was  always  less  than  0.1%. 

The  quantities  presented  are  the  azimuthal ly  averaged  mean  intensity 
(I)  and  vertical  flux  (F)  both  at  the  top  and  bottom  of  the  atmosphere. 

The  results  of  the  variational  technique  without  iterations  and  the  two- 
stream  method  (Liou,  1973)  are  also  included  in  the  tables  for  the  pur- 
pose of  comparison. 

Generally  speaking,  errors  in  the  variational  technique  lie  within 
2%  for  x*  = .25,  10%  for  t*  = 1 and  20%  for  x*  = 4.  On  the  other  hand, 
errors  in  the  two-stream  method  can  get  higher  than  50%  for  x*  - .25, 

30%  for  x*  = 1 and  20%  for  x*  = 4. 

From  these  results,  it  is  obvious  that  the  variational  technique 
is  much  more  accurate  than  the  widely  used  two- stream  method  for  the 
anisotropic  scattering  case  considered  (p(w,ii')  = u>q(1  + xwu')).  Further- 
more, once  the  algorithm  is  available,  the  computational  time  is  also 
minimal.  Therefore,  we  can  conclude  that  we  have  developed  a far  more 
powerful  technique  in  dealing  with  first  order  anisotropic  scattering. 

Inhomogeneity  is  an  important  aspect  of  the  radiative  transfer  prob- 
lem. In  most  approaches,  it  involves  matching  boundary  conditions  for 
every  region  of  different  ooo,  which  could  be  a tedious  job.  Inhomogen- 
eity can  be  directly  incorporated  in  the  variational  approach  through 
step  function  approximations  without  introducing  extra  computational 
complexities.  For  instance,  the  two-step  function  is  sufficient  for 
handling  a two-layer  cloud  problem.  The  capability  of  the  variational 
technique  to  represent  vertical  inhomogeneities  thus  provides  a great 
simplicity  in  treating  cloudy  atmospheres. 

Figures  1 and  2 provide  examples  for  inhomogeneous  atmospheres. 
Profiles  of  the  source  function,  mean  intensity  and  diffused  flux  are 
presented  for  the  two  different  cases  considered.  Both  cases  deal  with 


, t 'l  ' v ? Vk’A 

* 


"rlLf  '*£&**'  S :VV-7' 


TABLE  4 


CASES 

OF  ANISOTROPIC 

SCATTERING  WITH 

PHASE 

FUNCTION 

p(y.p')  = u>0(l  + 

xpvi')  and  t*  = 

.25,  yQ 

= 1. 

Exact 

Variational 

2-Strcam 

A% 

F+(0) 

0.63S2(-1) 

0. 6365 (- 1 ) 

0.2% 

0. 5538  (- 1) 

12.8% 

H 

r-* 

F+(t*) 

0.1573(0) 

0.1576(0) 

0.2% 

0. 1658(0) 

5.4% 

o 

II 

1(0) 

0. 5561 (-1) 

0. 5587  (-1 ) 

0.5% 

0. 2398  (- 1 ) 

55.9% 

3 

X 

T( T*) 

0.9573(-l) 

0.9025(01) 

0.5% 

0. 7180(-1) 

25.0% 

F+(0) 

0. 8818C-1) 

0. 881 7 ( - 1 ) 

0.0% 

0. 8426(-l) 

4.4% 

H 

to 

F+(t*) 

0.1327(0) 

0. 1 330  (- 1 ) 

0.2% 

0.1369(0) 

3.2% 

o 

n 

1(0) 

0.6680(-l) 

0. 6685(-l) 

0.1% 

0.  364 8 (-1) 

45.4% 

3 

B 

T(t‘) 

0. 8455 (- 1 ) 

0. 8528  (- ? ) 

0.9% 

0. 5930(-l) 

29.9% 

Ft(0) 

0.1549(0) 

0.1545(0) 

0.3% 

0.1612(0) 

4.1% 

H 

H 

1 

F+(t*) 

0.6602(-l) 

0.6674 (-1) 

1.1% 

0.6002  (-1 ) 

0.0% 

o 

II 

1(0) 

0.9706(-l) 

0. 9652 (- 1 ) 

0.6% 

0.6979(-l) 

28.1% 

I(T*) 

0.5428(-l) 

0. 5560(-l) 

2.4% 

0. 2599(-l ) 

52.1% 

to 

F + (0) 

0. 2288 (- 1 ) 

0. 2286(-l) 

0.1% 

0. 2293(-l) 

0.2% 

rH 

' F+(t*) 

0.6809(-l) 

0.6822(-l) 

0.2% 

0. 7603(- 1) 

11.7% 

II 

o 

M 

T(0) 

0.  2 1 79 (- 1 ) 

0-2172 (-1) 

0.3% 

0.9930 (-2) 

54.4% 

3 

X 

I(T*) 

0.4068(-l) 

0.4093(-l) 

0.6% 

0. 3292  (- 1 ) 

19.1% 

t/> 

FHO) 

0. 3457 (-1) 

0. 3464  (- 1 ) 

0.2% 

0. 3665(-l) 

6.0% 

to 

F+(t*) 

0. 5655(-l) 

0. 5640 (-1) 

0.3% 

0. 6238(-l) 

10.3% 

H 

o 

HO) 

0. 2689 (-1) 

0. 2702 (-1) 

0.5% 

0. 1587 (-1) 

41.0% 

3 

T(t*) 

0. 3579  (- 1 ) 

0. 3550(-l) 

0.8% 

0. 2701  (-1) 

24 . 5% 

to 

Ft(0) 

0.6810(-1) 

0. 6791  (-1) 

0.3% 

0. 7539  (- 1 ) 

10.7% 

H 

1 

F+(t*) 

0. 2315(-1) 

0.2337(-l) 

1.0% 

0. 2386(-l ) 

3.1% 

N 

O 

N 

T(0) 

0.4190(-l) 

0.4160(-1) 

0.7% 

0.3264(-l) 

22.1% 

3 

X 

Hi*) 

0.2075 (-2) 

0.2116C-1) 

1.9% 

0. 1033 (- 1 ) 

50.2% 

34 
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TABLE  5 


CASES  OF  ANISOTROPIC  SCATTERING  WITH  PHASE  FUNCTION 
p(p,p')  = u)q(1  + xpp')  and  t*  = 1,  = 1. 


r-H 

II 

o 


F (0) 

F (T*) 

T(0) 

T(x*) 


F (0) 

F (T.*) 
1(0) 
W) 


F (0) 

F (I*) 

T(0) 

1(7*) 


F (0) 

F (C*) 
1(0) 
T(-c*) 


Exact 


0.2326(0) 
0.  394 2(0) 
0.1413(0) 
0.  1944(0) 


0.2895(0) 
0.3374(0) 
0. 1661  (0) 
0. 1696(0) 


0.4181(0) 

0.2086(0) 

0.223(0) 

0.1133(0) 


0. 1404(0) 
0.4165(-1) 
0. 7375 (- 1 ) 
0.2495(-l) 


Variational 


0.2322(0) 
0. 3999(0) 
0.1428(0) 
0.2020(0) 


0.2869(0) 
0.3452(0) 
0. 1659(0) 
0.1789(0) 


0.4102(0) 

0.2220(0) 

0.2177(0) 

0.1270(0) 


A% 


0.2% 

1.4% 

0.4% 


1.9% 

6.4% 

2.1% 

12.1% 


0.1.361(0) 
0.4495(-l) 
0. 7023  (- 1 ) 
0. 2814 (-1) 


3.1% 

7.9% 

4.8% 

12.8% 


2-Stream 


0.2193(0) 
0.4128(0) 
0. 9498 (- 1 ) 
0.1787(0) 


0.2848(0) 

0.3473(0) 

0.1233(0) 

0.1504(0) 


0.4285(0) 
0.2036(0) 
0. 1856(0) 
0. 881 6 (-  1 ) 


F (0) 

0.  51 19 (- 1 ) 

0. 501 0 (- 1 ) 

2.1' 

i 0 . 5892 (- 1 ) 

f Cl*) 

0. 1205(0) 

0. 1236(0) 

2.6' 

6 0.1390(0) 

1(0) 

0.  364 3 (- 1 ) 

0. 351 6 ( - 1 ) 

3.5' 

i 0.2551 (-1 ) 

T(-t*) 

0.5721  (-1) 

0.6000(-l) 

4.9' 

6 0. 6021 (- 1) 

F (0) 

0.7591  (-1) 

0. 7402  (- 1 ) 

2.55 

6 0.8699 (- 1 ) 

F (X*) 

0.-9844  (-1) 

0.1016(0) 

3.  2! 

i 0.1138(0) 

T(0) 

0.4676(-l) 

0.4491 (-1) 

4.0? 

0.  3767  (-1 ) 

I(Tf) 

0.4817  (-1) 

0.5107 (-1) 

6. 0‘ 

0.4929  (- 1 ) 

0. 1S89 (0) 
0. 5026 (- 1 ) 
0. 6880 ( - 1 ) 
0.  21 76 (- 1 ) 


5.7% 
4.7 
32 . 8 
8.  1 


15.1% 

15.4% 

30.0% 

5.2% 


14.6% 

15.6% 

19.4% 

2.3% 


13.2% 

20.7% 

6.7% 

12.8% 


TABLE  6 


1 


CASES  OF  ANISOTROPIC 

SCATTERING  WITH  PHASE 

FUNCTION 

P(u,v')  = wo(l  + 

xpp')  and  t*  = 

4 and  y 

0 

= 1. 

Exact 

Variational 

A”; 

2-St  ream 

A% 

]••+«)) 

0.5728(0) 

0.5718(0) 

0.  2% 

0.5893(0) 

2.9% 

H 

F*(t*) 

0. 3831(0) 

0.4099(0) 

6.9°i 

0.3924(0) 

2.3% 

II 

ll 

T(0) 

0.2902(0) 

0.2988(0) 

3.0 % 

0.2552(0) 

12.1% 

U 

3 

X 

Ur*) 

0. 1674(0) 

0. 1897(0) 

13.3% 

0. 1699(0) 

1.5% 

F+(0) 

0.6286(0) 

0.6256(0) 

0.5% 

0.6503(0) 

3.5% 

H 

to 

F+(t*) 

0.3256(0) 

0.3561 (0) 

9.4% 

0.3314(0) 

1.8% 

II 

II 

T(0) 

0. 3143(0) 

0.3213(0) 

2.2% 

0.2816(0) 

10.4% 

3° 

X 

T(t*) 

0.  1424(0) 

0.1668(0) 

17.1% 

0.1435(0) 

0.8% 

F+(0) 

0. 7293(0) 

0.7224(0) 

0.9% 

0.7581(0) 

3.9% 

H 

r-4 

F+(T*) 

0.2201(0) 

0.2592(0) 

17.8% 

0.2236(0) 

1.6% 

II 

II 

T(0) 

0.3577(0) 

0.3617(0) 

1.1% 

0.3283(0) 

8.2% 

o 

3 

X 

T(t*) 

0.9681 (-1) 

0.1253(0) 

29.4% 

0.9683 (- 1 ) 

0.0% 

F+(0) 

0 . 61 22  (- 1 ) 

0.  5701  (- 1 ) 

6.9% 

0. 7348(- 1) 

20.0% 

»n 

F+U*) 

0.  24 1 7 (- 1 ) 

0. 2758(-l) 

14.1% 

0.2371 (-1) 

1.9% 

n 

e-^ 

II 

T(0) 

0. 4040  (- 1 ) 

0.  3626 (- 1 ) 

10.2% 

0. 3182  (- 1 ) 

21.2% 

o 

3 

X 

I(T*) 

0.96Q0(- 1) 

0. 1151  (-1) 

19.9% 

0. 1026  (-1 ) 

6.9% 

FUO) 

0.8927 (-1) 

0.8221 (-1) 

7.9% 

0.1041(0) 

16.6% 

to 

I/) 

F+(t*) 

0. 1839 (-1) 

0.2129(-1) 

15.8% 

0. 1832(-1) 

0.4% 

n 

II 

1(0) 

0. 5194  (- 1) 

0.4617  (-1) 

11.1% 

0.4508  (- 1 ) 

13.2% 

o 

3 

X 

T(I*) 

0. 7407 (-2) 

0.  9068(-2) 

22.4% 

0. 7931 (-2) 

7.1% 

Ft  (0) 

0.1590(0) 

0.1437(0) 

9:6% 

0.1792(0) 

12.7% 

»/> 

Ft  (r*) 

0.6655 (-2) 

0.  8427  (- 2) 

26.6% 

0. 7544 (-2) 

13.4% 

N 

1 

II 

T(0) 

0. 8076 (- 1) 

0. 7035 (-1) 

12.9% 

0.7761  (-1) 

3.9% 

3° 

X 

I(T*) 

0.2981 (-2) 

0.4094 (-2) 

37.3% 

0.3266 (-2) 

9.6% 

Figure  IB.  Mean  Intensity  (I)  and  Diffused 
Flux  (F)  Derived  from  Source 
Function  in  Figure  1A. 


an  atmosphere  of  total  optical  depth  (x*)  of  1 divided  into  two  regions 
such  that  from  x=0  to  0.5,  w = 1 and  from  x = 0.5  to  1.0,  w = .5.  In 
Figure  1,  X = 1 and  0.5  for  the  respective  regions  and  in  Figure  2, 

X = 0.5  and  1.  Figures  1A  and  2A  are  the  source  function  profiles. 

Figures  IB  and  2B  are  corresponding  profiles  of  the  mean  intensity  (I) 
and  diffused  flux  (F). 

Joseph  et  al.  (1976)  developed  a Delta-Eddington  approximation  to 
transform  the  Henyey-Greenstein  phase  function  into  the  form  of  1 + Xcos0 
(X  = 3g) . We  adopted  this  transformation  to  solve  for  different  atmo- 
spheric conditions  for  which  the  Henyey-Greenstein  phase  function  is  used. 
The  Henyey-Greenstein  phase  function  p(^  G(cos0)  is  defined  as 

PH-G  ^COS  °')  = (l+gz  - 2gScos  0)^3  ^1') 

in  which  g = <cos0>  is  the  asymmetry  factor.  Different  g values  represent 
variously  peaked  phase  functions.  It  can  vary  from  -1  (complete  back 
scattering)  through  0 (isotropic  scattering)  to  g = 1 (fully  foward  scat- 
tering). The  Henyey-Greenstein  phase  function  is  a useful  phase  function 
because  it  applies  to  a variety  of  different  atmospheres.  Haze  has 
asymmetry  factor  g ^ 0.75.  Cloud  droplets  and  aerosol  particulates  have 
g factors  between  0.85  and  0.95  for  most  of  the  solar  spectrum. 

The  Legendre  polynomial  expansion  of  the  H-G  phase  function  is  shown 
to  have  the  simple  form  (van  de  Hulst,  1968): 

CO 

PH  G(cos  0)  ~ l (21  + l)gzPjl(cos  0)  (52) 

1 = 0 

where  P^fcosO)  is  the  Legendre  polynomial  expansion.  The  Delta-Eddington 
approximation  is  a special  combination  of  the  Eddington  and  forward-peak- 
truncation  approximations.  It  can  be  written  as 

PH-G(C0S  0)  * P6-edd  (cos  0) 

= 2f  6 (1-cos  0)  + (1-f ) (l+3g 'cos  0)  (53) 

OO 

= 1 + 3g  cos  0+1  (21  + l)g2Pjl(cos  0) 

1 = 0 

where  f is  the  fractional  scattering  into  the  forward  peak  and  g'  is  the 
asymmetry  factor  of  the  truncated  phase  function,  f and  g are  shown  to  be 


and 


(54) 

(55) 


f = 


2 

g 


g' 


= _S_ 
1+g 


The  Delta-Eddington  phase  function  agrees  with  that  of  the  H-G  out 
to  three  terms  and  their  difference  is 

oo 


P6-edd(cOS  0)  - PH-G(cos  0) 


l (2t-l)(g2-gV»(cos  0) 
4=3 


(56) 


The  error  in  the  phase  function  tends  to  grow  smaller  as  g •+■  1 . At  g = 1 
(total  forward  scattering),  the  Delta-Eddington  approximation  is  actually 
exact.  By  redefining 


<(t) 


(l-f)M0(T) 

l-fw0(x) 


(57) 


T'  = (l-f«00(T))T  (58) 

and  I'(t'.p)  = I(1~ p^0(V) ,lJ)  = (59) 

we  obtain  the  transformed  set  of  the  radiative  transfer  equation  and  the 
coupled  integral  equations  for  the  source  function. 

Tables  7 and  8 show  results  of  some  of  the  tested  cases.  The  doub- 
ling method  was  assumed  to  yield  exact  solution  for  comparison  on  the 
accuracies.  Asymmetry  factor  g = 0.75  is  used.  u>o  = 1 and  0.8  represent 
fully  and  partially  scattering  atmospheres,  t*  = 0.25,  1 and  4 are 
typical  of  thin,  medium,  and  thick  atmospheres.  = 0.9  and  0.5  are 
chosen  for  near-normal  and  grazing  solar  incidences.  The  quantities 
presented  are  diffused  fluxes  at  the  top  and  bottom  of  the  atmosphere 
(F+ (0)  and  F+(i*)  respectively). 

The  results  show  that,  despite  the  truncation  in  the  phase  function 
the  percentage  error  is  relatively  small.  It  is  thus  conceivable  to 
employ  our  method  in  real  cloud  conditions  where  anisotropy  is  non-negli- 
gible.  In  the  next  section  we  will  discuss  some  examples  of  cirrus  clouds 
in  different  wavelength  ranges. 
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TABLE  7 


T* 

% 

P 

Exact 

Variational 
- iterative 

Error  % 

0.9 

P+(0) 

R(t*) 

0 . 2025  (-  1 ) 
0.1980 

0 . 2375  (- 1 ) 

0. 1945 

1 7 . .3* 
1.8°. 

U • 4 J 

0.5 

F+(0) 

0. 3590  (- 1 ) 
0.1608 

0.  .3455  (- 1) 
0.1620 

3.4% 

0 . 7% 

1 

0.9 

Ft  (0) 

F+  (t*0 

0 . 8703  (-  1) 
0.5167 

0.9120(-1) 

0.5076 

8.2". 

1.8% 

0.5 

Ft(0) 

F+(t*) 

0.1202 

0.5121 

0.1135 

0.3189 

S . 6% 
2.1% 

4 

0.9 

Fi(0) 

F-t  (i*) 

0.3134 

0.5760 

0.3145 

0.5673 

0 . 3% 
1.5% 

0.5 

Ft(0) 

FM'*) 

0.2597 

0.2402 

0.2546 

0.2421 

1.9% 

0.8% 

Accuracy  of  the  VI  Technique  (Doubling  Method  is  used 
to  obtain  exact  solutions.  Phase  function  is  that  of 
the  H-G  with  Delta-Eddington  transformation  and  o>rt  = 

, g = 0.75.) 


3.3.2  The  Effects  of  Multiple  Scattering  on  Microwave  Brightness 
Temperature  Observations 

Three  different  tropical  cloud  (with  rain)  models  are  considered  in 
the  simulation  of  microwave  brightness  temperature  values,  one  being 
single  layer  and  the  other  two  multilayer.  An  ocean  background  with  a 
surface  temperature  of  300°K  is  assumed. * Surface  reflectivities  of  0.6 
for  19.35  GHz  and  0.552  for  37  GHz  are  used. 

The  optical  depth  and  single  scattering  albedo  values  were  calcu- 
lated from  the  given  profiles  of  volume  extinction  and  absorption  coef- 
ficients of  water  vapor,  oxygen  and  cloud  droplets.  The  optical  depth 
and  single  scattering  albedo  profiles  obtained  are  then  approximated  by 
different  step  functions  (average  step  size  is  At  'v  0.2)  where  each  step 
is  assigned  a fixed  value  of  single  scattering  albedo.  The  source  func- 
tion profile  is  then  obtained  by  applying  the  variational-iterative  tech- 
nique. The  brightness  temperature  (proportional  to  the  outgoing  radiance) 
is  thus  calculated  from  (16)  for  different  look  angles.  Tables  9 to  11 
show  the  resultant  brightness  temperature  for  each  of  the  cloud  cover 
cases  considered. 

The  contribution  of  emission  from  the  atmosphere  and  clouds  (first 
term  in  (16))  is  always  the  least  for  y = 1 (normal  look  down  angle),  and 
increases  for  decreasing  y.  The  opposite  is  true  for  surface  emission 
and  reflection  (second  and  third  terms  in  (16)).  Furthermore,  the  contri- 
bution from  the  background  surface,  both  emission  and  reflection,  drops 
exponentially  as  a function  of  the  total  optical  thickness.  As  a result, 
for  an  optically  thick  atmosphere  (e.g.  Table  11B)  the  total  brightness 
temperature  is  the  greatest  at  y = 1 because  surface  emission  and  reflec- 
tion terms  are  appreciable  only  near  y * 1 and  angular  dependence  on  the 
atmospheric  emission  is  less  significant.  For  thinner  atmospheres  on  the 
other  hand,  the  total  brightness  temperature  is  the  least  at  y = 1 due 
to  the  relatively  small  contribution  of  the  atmospheric  emission  near 
y = 1 (e.g.  Tables  9,  10  and  11A). 

The  results  are  also  compared  with  those  from  the  ERT  Radiative 
Absorption  and  Scattering  Program  (RASP)  with  the  same  cloud  conditions 
applied.  The  main  difference  between  the  two  models  is  that  RASP  does 
not  handle  multiple  scattering  through  the  source  function;  the  source 
function  is  simply  the  Planck's  function  (J(t)  = B(T(t)))  in  which 
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TABLE  9 

One  layer  tropical  cloud  (0-4000  m) 

Liquid  water  content  = 0.045  gm/m3 
Mode  radius  = 400  v 
Surface  temperature  = 300° K 

A.  Frequency  19.35  GHz 
Surface  albedo  = 0.6 

Total  optical  depth  = 0.3  (0.14  without  cloud) 

Single  scattering  albedo  up  to  4000  m = 0.15  to  0.22 


Brightness  Temperature  Calculated: 


V 

Atm.  Emi . 

Surf.  Kef. 

Surf.  Emi . 

Total 

0.2 

201.9 

14.2 

26 . 8 

24  2.9 

0.4 

L3S.  3 

30.  2 

56.7 

225.  1 

0.6 

103.4 

58 . 7 

72. 8 

21.4.9 

o.s 

82.3 

4 3.9 

82.5 

208.7 

1 .0 

68.  3 

47. 3 

SS.9 

201 . 5 

Unit:  °K 


B.  Frequency  37.0  GHz 

Surface  albedo  = 0.552 

Total  optical  depth  = 0.5  (0.16  without  cloud) 

Single  scattering  albedo  up  to  4000  m * 0.37  to  0.46 


Brightness  Temperature  Calculated: 


P 

Atm.  Emi . 

Surf.  Kef. 

Surf. 

Em  i . 

Total 

0.  2 

205.8 

6.0 

11. 

0 

222.8 

0.4 

163.6 

21.1 

38 

5 

223.2 

0.6 

130.7 

32.0 

58 

4 

221.1 

0.8 

107.9 

39.4 

71 

4 

219.2 

Lljl_ 

91.6 

44.6 

81 

5 

217.7 

Unit:  #K 


1 

( 


46 


T 


TABLE  10 


SBZS  MOB  1$  BBST  ^UALmrPB^jp^m 
UtOM  COPY  TURHISUD  lOPM  — 


Five  layer  tropical  cloud  I 

(0-2000  m,  2000  m - 4000  m,  4000  m - 6000  m,  6000  - 8000  m, 

8000  - 10,000  m) 

Liquid  water  content  = 0.060,  0.050,  0.030,  0.020  and  0.020  gm/m 
respectively 

Mode  radius  = 450,  400,  350,  300,  200  p respectively 
Surface  Temperature  = 300°K 

A.  Frequency  19.35  GHz 
Surface  albedo  = 0.6 

Total  optical  depth  = 0.24  (0.15  without  cloud) 

Single  scattering  up  to  10,000  m = 0.009  to  0.044 


Brightness  Temperature  Calculated: 


p 

Atm.  Emi . 

Surf.  Kef . 

Surf.  Emi . 

Total 

0.  2 

192.  r 

17.3 

36. 1 

245.5 

0.4 

174.9 

31.4 

65.9 

222. 2 

0.6 

91.5 

38. 4 

80.4 

210.  5 

0.8 

72.0 

42.5 

88.9 

203.4 

1 . 0 

59.  5 

45.  1 

94.4 

198.  8 

Unit:  *K 


B.  Frequency  37.0  GHz 

Surface  albedo  = 0.552 

Total  optical  depth  = 0.52  (0.17  without  cloud) 

Single  scattering  albedo  up  to  10,000  m = 0.03  to  0.23 


Brightness  Temperature  Calculated: 


P 

Atm.  Emi . 

Surf.  Kef. 

Surf.  Emi . 

Total 

0.2 

242.8 

6.3 

10.0 

259.1 

0.4 

193. 1 

23.  2 

36.5 

252.9 

0.6 

154.5 

35.8 

56.  5 

246.8 

0.8 

127.7 

44.5 

70.  2 

242.3 

1.0 

108.4 

50.7 

79.  9 

239.0 

Unit:  °K 
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TABLE  11 


Five  layer  tropical  cloud  II  j 

(0  - 2000  m,  2000  m - 4000  m,  4000  m - 6000  m,  6000  m - 8000  m, 

8000  - 10,000  m)  3 j 

Liquid  water  content  = 0.360,  0.300,  0.180,  0.130,  0.100  gm/m, 

respectively  . I 


Mode  Radius  = 450,  400,  350,  300,  200  y,  respectively 
Surface  temperature  = 300° K 

A.  Frequency  19.35  GHz 
Surface  albedo  = 0,6 

Total  optical  depth  = 0.71  (0.15  without  cloud) 

Single  scattering  albedo  up  to  10,000  m = 0.01  to  0.086 


Brightness  Temperature  Calculated: 


Surface  albedo  = 0.552 

Total  optical  depth  = 2.24  (0.17  without  cloud) 

Single  scattering  albedo  up  to  10,000  m = 0.03  to  0.32 


Brightness  Temperature  Calculated: 


f 


V 


I 

I 


scattering  cross  sections  are  included  for  calculation  of  the  optical 
depth.  Another  difference  is  the  handling  of  the  surface  reflection 
terms.  In  the  VI  program,  the  Lambertian  surface  approach  is  adopted. 
Reflectivity  of  the  surface  is  defined  as  the  ratio  between  the  incoming 
and  outgoing  fluxes  at  the  surface  and  the  reflection  is  isotropic  (i.e. 
no  angle  preference  for  the  outgoing  radiance)  (equ.  (14)).  In  RASP,  the 
specular  reflection  condition  is  used,  or 

I(t*,u)  = R I(t*,-u)  (60) 

As  a result,  because  the  atmospheric  emission  term  is  generally  higher 
for  smaller  y's,  at  u = 1 (normal  look  down  angle)  the  Lambertian  model, 
which  takes  incoming  radiation  from  all  angles  into  consideration, 
gives  a higher  surface  reflection  term  than  that  of  the  specular  ref*  . c- 
tion. 

Comparisons  are  carried  out  for  the  same  clouds  at  p = 1,  the  angle 
which  RASP  operates  specifically.  Specular  reflection  at  the  surface 
for  the  VI  model  is  included  to  demonstrate  the  difference  between  dif- 
ferent reflection  models.  The  results  of  different  approaches  are  shown 
in  Table  12. 

It  is  shown  from  the  last  two  columns  in  Table  12  that  the  differ- 
ence in  the  brightness  temperatures  between  specular  and  Lambertian  sur- 
face reflections  generally  lies  between  15  and  20  degrees  except  for  the 
thick  atmosphere  case  (t*  = 2.24)  where  the  difference  is  only  4 degrees 
because  in  this  case  the  angular  dependence  on  the  atmospheric  emission 
term  is  not  significant  (Table  11B). 

Results  of  RASP  and  VI  techniques  are  also  compared  for  the  same 
surface  model.  It  is  shown  (Table  12)  that  the  two  models  agree  quite 
well  for  small  single  scattering  albedos  (u)q  _<  0.05).  When  appreciable 
scattering  is  present,  RASP  tends  to  overestimate  the  brightness  temper- 
ature for  thinner  atmospheres  because  the  scattering  term  (second  term 
in  (15))  is  less  than  u>  B.  On  the  other  hand,  RASP  also  may  underesti- 
mate the  brightness  temperature  when  the  optical  depth  is  higher  for  the 
more  complex  conditions  specified  in  Table  11. 
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TABLE  12 


Summary  of 

Brightness  Temperature 

no- 

Atmos 

pheric  Conditions 

VI 

VI 

T* 

<i>0 

R 

■BE*/  . • ■ H 

(specular  ref) 

(Lambertian) 

0.3 

0.15-0.22 

195.4 

186.1 

204.5 

0.5 

0.31-0.46 

0.552 

232.4 

198.3 

217.7 

0.24 

B 

182.2 

181.2 

198.8 

0.52 

1 ' . 

209.2 

221.5 

239.0 

0.71 

B 

232.9 

237.9 

254.3 

2.24 

1 

0.552 

233.5 

264.1 

268.9 

V = 1 


C 


5 


3.3.3  Anisotropic  Scattering  through  Cirrus  Cloud  Layers  in  the 
Visible  and  Infrared  Regions 

The  major  difficulty  of  scattering  computation  for  cirrus  clouds 
arise  from  the  random  distribution  of  the  highly  nonspherical  ice  cry- 
stals in  the  cloud  layers.  Little  information,  either  theoretical  or 
observational,  on  the  composition  of  the  cirrus  clouds  is  known.  However, 
though  geometrically  thin  (<500  m) , they  are  usually  optically  thick 
which  makes  the  determination  of  the  atmospheric  structure  from  orbiting 
meteorological  satellites  difficult. 

A preliminary  study  of  the  scattering  problem  through  cirrus  clouds 
has  been  performed.  We  adopted  Liou's  (1972)  computations  of  the  phase 
function,  the  single  scattering  albedo  and  the  extinction  cross-section 
(Bext)  of  cirrus  cloud  layers  for  some  sample  wavelengths  (0.7  p,  3 p, 
and  10  p) . The  assumption  is  that  the  ice  crystals  in  cirrus  clouds  can 
be  approximated  by  long  circular  cylinders  which  are  randomly  oriented 
in  space.  Table  13  is  a summary  of  the  optical  properties  (ujq , g,  8 t) 
of  randomly  oriented  ice  crystals  at  the  sample  wavelengths.  We  then 
chose  three  cirrus  cloud  layers  with  different  thicknesses:  100  m for  a 
very  thin  layer,  200  m for  a medium  thin  layer,  and  500  m for  a relatively 
thick  cirrus  cloud  layer.  The  total  optical  depth  of  each  layer  is  simply 
'the  product  of  the  volume  extinction  cross-section  (8  t)  and  the  geometric 
thickness . 

Tables  14,  15,  and  16  are  results  for  A = 0.7  p,  3 p and  10  p respec- 
tively. Cosines  of  solar  zenith  angles  (pQ)  of  0.9  and  0.5  are  chosen 
as  cases  for  near-normal  and  grazing  solar  incidences.  Quantities  pre- 
sented are  reflection  (r)  and  transmission  (t)  as  defined  as 

r = itF+(o)/ttpoFo  (61) 

and  t = ttFI  (t*)/ttpoFo  ♦ exp(-T*/pQ)  (62) 

where  irF  is  the  unattenuated  solar  flux.  Results  with  isotropic  scat- 
tering (g  = 0)  for  the  same  cases  are  also  included  to  demonstrate  the 
effect  of  forward  scattering  of  the  ice  crystals. 

As  expected,  g > 0 (forward  scattering)  always  increases  the  trans- 
mission and  decreases  the  reflection  as  compared  to  g = 0.  It  is  also 
shown  that  in  the  visible  range,  most  of  the  light  gets  transmitted  through 


51 


TABLE  15 

l 


Thickness 

Uo 

r 

r1 (isotropic) 

t 

t* (isotropic) 

100  m 

0.9 

0.3359(-l) 

0.1074 

0.2386 

0.1706 

0.5 

0.7217  (-1) 

0.1310 

0.9527 (-1) 

0.7416(-1) 

200  m 

0.9 

0.3415(-1) 

0.1096 

0 . 6760  C - 1 ) 

0.2623(-l) 

0.5 

0. 7354 (-1) 

0.1320 

0.2046(-l) 

0,8386(-l) 

500  m 

0.9 

0.3426(-l) 

0.8783(-l) 

0. 5144 (-3) 

0 .3317 (-3) 

0.5 

0.7441  (-1) 

0.9490(-l) 

0 . 1081 (-3) 

0.8526(-4) 

Reflection  (r)  and  Transmission  (t)  through 
Different  Cirrus  Cloud  Layers  at  X = 3 u 


TABLE  16 


l 


Thickness 

% 

r 

r*  (isotropic) 

t 

t* (isotropic) 

100  m 

0.9 

0.9571 (-2) 

0.1074 

0.3052 

0.1706 

0.5 

0.2948(-l) 

0.1310 

0.1266 

0.7416(-1) 

200  m 

0.9 

0.9832C-2) 

0.1096 

0.8085(-l) 

0.2623(-l) 

0.5 

0.3030(-l) 

0.1320 

0. 1691 (-1) 

0. 8386 ( - 2) 

500  m 

0.9 

0. 9842(-2) 

0.8783(-l) 

0. 2384 (-2) 

0.3317 (-3) 

0.5 

0 . 3053 (- 1 ) 

0.9490(-l) 

0.2753(-3) 

0.8526(4) 

Reflection  (r)  and  Transmission  (t)  through 
Different  Cirrus  Cloud  Layers  at  X = 10  y 
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the  cloud  layer.  Furthermore,  since  w = 1,  a case  of  total  scattering, 
there  is  no  absorption  through  the  cloud  layer.  (The  fraction  of  absorp- 
tion a = 1 - r - t.)  However,  into  the  IR  region,  except  for  very  thin 
layers,  most  of  the  light  is  absorbed.  Reflection  or  transmission  is  at 
most  a few  percent. 
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4.  INFRARED  RADIANCE  SIMULATIONS  WITH  CLOUDS 


The  infrared  radiance  simulation  package  obtained  from  AFGWC  and 
modified  for  use  on  the  AFGL  CDC  6600  computer  was  further  modified  to 
incorporate  cloud  scattering  and  the  use  of  the  VI  technique  for  calcu- 
lating the  radiance  values.  The  program  was  modified  only  to  incorporate 
isotropic  scattering,  the  effects  of  anisotropy  being  of  second  order 
relative  to  the  effect  of  including  scattering  in  the  radiance  calcula- 
tions . 

4.1  Software  Modification 

The  simulation  program  RAPID  GABTAWF  (Weichel,  1976)  was  modified 
to  include  the  VI  algorithm  for  the  multiple  scattering  effect  rather 
than  the  numerical  quadrature  procedure  assuming  no  scattering  used  by 

t 

AFGWC.  The  transmittance  values  for  the  70  levels  of  the  AFGWC  integra- 
tion program  were  multiplied  by  tlje  additional  cloud  contributions  and 
the  resultant  weighting  function  profile  was  used  to  generate  a new  set 
of  heights  for  the  VI  subroutine.  The  profile  was  divided  into  roughly 
20  levels  of  nearly  equal  optical  depth,  each  level  being  approximated 
by  an  averaged  single  scattering  albedo.  These  levels  were  used  for  the 
multiple  scattering  computations  and  further  iterations  were  applied  for 
better  accuracies. 

4.2  Cloud  Simulations 

The  variational-iterative  technique  was  fully  tested  with  realistic 
atmospheres  and  clouds  using  the  CDC  6600  computer  facilities  at  AFGL. 

The  effects  of  clouds  were  simulated  using  two  cloud  models:  high- 
thin  cirrus  clouds  and  low  level  stratus  clouds.  These  cases  were  chosen 
because  they  represent  the  conditions  under  which  the  normal  cloud  screen- 
ing algorithms  would  fail.  Table  17  shows  the  properties  of  these  clouds 
including  the  cloud  type,  base  and  top  altitudes,  density  and  mode  radius 
of  the  cloud  droplets.  In  all  cases,  the  clouds  were  assumed  to  com- 
pletely fill  the  field  of  view.  A radiosonde  profile  from  Washington, 

D.C.  (spring)  with  surface  temperature  of  288°K  and  emissivity  of  one  was 
used.  The  refractive  indices  of  ice  and  water  were  provided  by  the  con- 
tract monitor,  Mr.  V.J.  Falcone. 
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TABLE  17 


Calculations  for  radiances  were  done  at  8 IR  channels  with  wave- 
lengths of  15,  14.8,  14.4,  14.2,  13.8,  13.4,  12,  and  18,7  microns, 
respectively.  The  7th  channel  (12  micron)  is  a window  channel  which  has 
high  transmittance  under  clear  atmospheric  conditions.  The  optical  pro- 
perties of  the  clouds  at  this  channel,  t ^ (total  optical  depth)  and  w 
(single  scattering  albedo)  are  also  presented  in  Table  17. 

Results  for  various  clouds  are  shown  in  Figures  3-8  for  all  eight 
channels.  Radiances  (B)  are  also  converted  to  the  equivalent  black  body 
temperature  (Tg)  defined  by  the  Planck  radiance  function: 

B = 2hcV5/(exp(hc/kATD)-l)  (63) 

D 

where  h is  Planck's  constant,  c is  the  speed  of  light  and  A is  the  wave- 
length. Computations  were  done  both  for  pure  absorption  (the  AFGWC 
model)  and  multiple  scattering  (VI). 

It  has  been  demonstrated  that  the  VI  technique  is  generally  fast  and 
relatively  accurate.  Tests  were  carried  out  using  different  layer  thick- 
nesses. It  was  shown  that  VI  deals  with  thin  atmospheres  the  best  (in 
contrast  to  the  Doubling  Method  and  Discrete  Ordinate  Method) . For 
thicker  atmospheres,  a larger  number  of  layers  and  more  iterations  were 
required  for  the  desired  accuracy.  It  is  also  noted  that  for  a discon- 
tinuous atmosphere  (for  example,  a scattering  cloud  imbedded  in  a large, 
purely  absorbing  atmosphere)  the  model  should  include  the  multiple  scat- 
tering (for  the  cloud)  and  pure  absorption  (for  the  rest  of  the  atmo- 
sphere) computations  simultaneously. 

Results  show  that  for  channels  1 through  5 the  atmosphere  above  the 
cloud  is  sufficiently  opaque  to  mask  the  effects  of  multiple  scattering 
within  the  cloud.  In  other  words,  a pure  absorption  treatment  is  ade- 
quate and  the  general  assumption  of  brightness  temperature  equal  to  the 
cloud  top  temperature  as  used  in  Section  2 is  valid  for  these  channels. 

Secondly,  for  all  low  stratus  clouds,  no  scattering  treatment  is 
necessary  except  in  the  window  channel  (No.  7).  This  is  verified  by  the 
facts  that  in  these  cases  (1)  there  is  little  difference  in  radiance 
values  between  pure  absorption  and  scattering  computations,  and  (2)  the 
equivalent  black  body  temperature  is  always  below  the  cloud  top  temper- 
ature indicating  the  already  low  transmittance  values  before  entering 
the  cloud  top. 
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Brightness  Temperature 
( indicates  the 


Brightness  Temperatures 
( indicates  the  cL 


Cirrus  clouds,  on  the  other  hand,  need  multiple  scattering  treatment 
due  to  their  high  altitudes  and  the  more  transparent  optical  properties 
(lower  optical  depth  values).  It  is  shown  that  up  to  a 5%  difference  in 
radiances  (or  <5°  cooling  in  black  body  temperature)  can  result  not  only 
in  channel  no.  7 but  also  in  channels  no.  6 and  8. 

Another  aspect  of  cirrus  clouds  should  be  mentioned.  Due  to  the 
existence  of  ice  crystals,  anisotropic  scattering  should  be  taken  into 
account.  Treatment  of  forward  scattering  properties  of  the  ice  crystals 
need  to  be  included.  In  this  section  we  classified  the  cases  of  cloud 
models  in  the  IR  region  that  need  the  inclusion  of  scattering  effects. 
Further  effort  should  emphasize  the  need  for  an  anisotropic  scattering 
treatment  for  cirrus  clouds  with  the  use  of  the  complete  phase  functions 
for  ice  crystals.  Such  investigation  would  improve  the  accuracy  of  pre- 
sent models  for  thin  cirrus  problems. 
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5.  PARTIAL  CLOUDINESS  IN  THE  FIELD-OF-VIEW 

The  variational-iterative  procedure  for  obtaining  a solution  to  the 
radiative  transfer  equation  was  developed  for  application  to  problems 
with  axial  symmetry  about  the  zenith  direction  in  a plane-parallel  scat- 
tering atmosphere.  This  limitation  is  inherent  in  most  of  the  rapid 
techniques  for  obtaining  a solution  to  the  radiative  transfer  equation 
In  the  case  of  the  variational-iterative  technique,  the  restrictions  may 
be  relaxed  but  only  at  the  expense  of  significantly  increasing  the  dimen- 
sions of  the  matrices  that  must  be  inverted  to  obtain  a solution.  In 
special  cases,  this  may  provide  more  insight  into  the  nature  of  multiple 
scattering  problems  but,  in  the  general  case  of  simulating  brightness 
fields  for  horizontally  homogeneous  clouds,  the  application  of  this  or 
any  of  the  other  plane-parallel  atmosphere  procedures  is  impractical. 

The  problem  of  obtaining  temperature  profiles  from  brightness  mea- 
surements in  the  presence  of  clouds  was  reviewed  in  Section  2.3.  The 
procedures  used  were  based  upon  the  assumption  that  the  region  with 
clouds  could  be  subdivided  into  subregions  with  and  without  clouds  and 
that  the  radiative  transfer  problem  could  be  independently  solved  for 
each  subregion.  The  total  radiance  observed  at  the  satellite  then  could 
be  obtained  by  summing  the  contributions  from  each  region.  Since  the 
different  regions  contributed  in  proportion  to  their  areas,  the  area 
weighted  results,  equation  (4),  is  obtained. 

The  response  of  a satellite-borne  radiometer  to  a horizontally 
inhomogeneous  scene  (surface  brightness  field)  viewed  through  a horizon- 
tally inhomogeneous  non-scattering  atmosphere  may  be  represented  by  the 
generalization  of  equation  (1): 

P 

* s 

V //s^n){BJT(ps,xs,Ys)]0v(Ps.VV  - / WP.x.y)]  • 

surface  0 (64) 

■d-G-vl£— ’ y-“  dp>  dX  dY 
dP  * s s 

where  5,n  represent  the  angle  subtended  by  the  point  X ,Y  as 
viewed  from  the  satellite, 

g(£,n)  is  the  intensity  response  of  antenna  gain  of  the  opti- 
cal or  microwave  system  in  the  direction  £,Y 


i 

I 


I 


i 
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and  x,y  represent  the  horizontal  position  of  the  ray  between  the 

satellite  and  P , Xg,  Yg  at  the  pressure  level  P. 

If  the  intensity  response  function,  g(£;,n)  is  a delta  function,  then  only 
radiation  from  a single  direction  So,no  is  observed  and  equation  (1) 
obtains.  If  the  response  function  is  narrow  so  the  changes  in  the  inte- 
grand  specified  by  variations  in  Xs>  Y , X,  and  Y are  small  in  compari- 
son with  changes  in  g(n,£),  equation  (1)  still  obtains.  However,  if  the 
response  function  has  a large  field-of-view  with  a uniform  response  in 
each  direction,  g(£,o)  = constant,  then  equation  (64)  corresponds  to  the 
area  weighted  response,  equation  (4) . 

The  use  of  equation  (64)  to  represent  the  area  weighted  cloud  con- 
tribution requires  the  additional  assumption  that  the  clouds  may  be 
locally  treated  as  plane-parallel  scattering  regions  and  the  relative 

contributions  of  the  edges  are  small  (locally  horizontally  homogeneous). 

v 

While  this  assumption  can  only  be  checked  using  the  full  radiative  trans- 
fer equation  and  horizontally  inhomogeneous  model  atmospheres,  the  uncer- 
tainties in  the  model  clouds  used  to  represent  the  atmosphere  may  be 
larger  than  the  effects  of  the  inhomogeneities.  In  the  end,  the  use  of 
equation  (64)  to  investigate  the  effects  of  clouds  depends  upon  assump- 
tions about  the  clouds  within  the  field-of-view  and  model  clouds  that 
are  relatively  much  larger  in  the  horizontal  than  vertical  may  be  postu- 
lated that  reasonably  fulfill  the  requirements  for  local  homogeneity. 
These  clouds  must  be  treated  as  effective  clouds  with  properties  inferred 
to  correspond  to  observations. 

• 


6.  CONCLUSIONS 


The  variational-iterative  (V-I)  procedure  for  the  solution  of  the 
radiative  transfer  equation  for  a plane-parallel  scattering  atmosphere 
(multiple  scattering)  was  modified  and  extended  for  use  in  problems 
involving  anisotropic  scattering.  The  V-I  technique  was  incorporated 
in  the  AFGL  simulation  package  for  the  DMSP  infrared  sounder.  Applica- 
tion of  the  technique  to  conditions  of  high  thin  cirrus  clouds  or  to  low 
stratus  clouds  showed  that  the  multiple  scattering  calculations  are 
required  to  adequately  simulate  cloud  effects  in  the  window  channel;  the 
effect  of  multiple  scattering  was  to  produce  a five  percent  difference 
in  the  simulated  radiance  value,  a value  large  compared  with  the  expected 
measurement  noise. 

A statistical  analysis  of  the  relative  merits  of  the  infrared  only, 
microwave  only,  and  infrared  plus  microwave  sounder  systems  revealed 
that  for  cloud  free  skies,  one  system  did  as  well  as  the  other.  In  the 
presence  of  clouds  and  for  microwave  and  infrared  sounder  systems  with 
identical  f ields-of-view,  the  microwave  sounder  system  was  superior  to 
the  infrared  system;  the  addition  of  infrared  data  to  the  microwave 
sounder  data  provided  new  information  for  the  assessment  of  cloud  cover 
and  height  but  no  new  information  for  the  deduction  of  temperature  pro- 
files. For  the  current  and  planned  microwave  sounder  systems,  the  much 
larger  field-of-view  of  the  microwave  system  can  be  a disadvantage. 

Using  current  procedures  for  the  evaluation  of  temperature  profiles  in 
the  presence  of  clouds,  the  profile  inferred  from  the  infrared  sounder 
is  assumed  to  apply  to  a horizontal  region  significantly  larger  than  the 
field-of-view  of  the  infrared  sensor.  When  the  microwave  sensor  field- 
of-view  matches  the  larger  horizontal  region  of  the  infrared  system  in 
the  presence  of  clouds,  the  microwave  system  is  superior  for  temperature 

profiling  and  the  infrared  system  should  be  used  in  combination  with  the 

* 

microwave  system  for  estimating  cloud  height  and  coverage. 
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